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ABSTRACT 
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Surface  velocity  fields  from  two  configurations 
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Los  Alamos  National  Laboratory  (LANL)  Parallel  Ocean  Program 
(POP)  model  are  compared  to  surface  velocity  data  from 
satellite-tracked  buoys  in  the  North  Atlantic.  Separate 
analyses  are  conducted  for  each  model  configuration.  In  the 
first  analysis,  output  from  a  1/6-degree,  20-level  model 
version  is  compared  with  five  years  (1993-1997)  of  drifter 
data,  based  on  both  Eulerian  and  Lagrangian  statistics.  In 
the  second  analysis,  newly-available  output  from  a  1/10- 
degree,  40  level  version  is  compared  to  a  two-year  subset 
(1993-1994)  of  the  data,  and  to  1/6-degree  output  over  the 
same  time  frame.  The  latter  comparison  is  based  on  Eulerian 
statistics  alone. 

The  five-year  comparison  shows  that  the  1/6-degree 
model  produces  inaccuracies  in  some  features,  and  generally 
underestimates  velocity  variance.  Modeled  Lagrangian  time 
scales  are  too  long,  while  the  length  scales  are  too  short. 
The  two-year  comparison  shows  that  at  the  higher  vertical 
and  horizontal  resolution  of  the  1/10-degree  model,  there  is 
a  striking  improvement  in  the  spatial  distribution  of  energy 
and  resolution  of  the  variance  field. 
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EXECUTIVE  SUMMARY 

To  meet  future  Navy  operational  ocean  forecasting 
requirements,  a  high-resolution,  coupled  air/ocean/ice  model 
with  data  assimilation  capability  will  be  needed.  Current- 
generation  ocean  general  circulation  models  have  evolved  to 
a  point  that  suggests  evaluating  them  for  this  role.  The  Los 
Alamos  National  Laboratory's  (LAND  Parallel  Ocean  Program 
(POP)  is  one  of  the  ocean  models  currently  under 
consideration  by  the  Office  of  Naval  Research  (ONR) .  The 
purpose  of  this  study  is  to  evaluate  the  ability  of  the  POP 
to  reproduce  surface  features  and  processes  of  interest  to 
the  Navy,  by  statistical  comparison  of  its  velocity  fields 
to  those  of  surface  drifters  in  the  North  Atlantic  Ocean. 

We  compare  output  from  a  1/6-degree,  20-level  version 
of  the  model  to  five  years  of  drifter  data  spanning  the 
period  January  1,  1993  to  December  31,  1997,  on  the  basis  of 
Eulerian  and  Lagrangian  statistics.  Trajectories  for  the 
Lagrangian  analysis  are  generated  by  numerical  integration, 
using  bilinear  interpolation  over  the  model's  discrete 
velocity  fields. 

To  assess  the  effect  of  increased  horizontal  and 
vertical  resolution  in  the  model,  we  use  Eulerian  statistics 
to  compare  two  years  of  newly-available  output  from  a 
1/10-degree,  40-level  model  run  to  the  first  two  years  of 
the  drifter  data. 
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The  Eulerian  analysis  of  the  1/6-degree  model  indicates 
that  this  version  performs  best  at  higher  latitudes  (roughly 
north  of  6  0°N) ;  the  northern  boundary  currents  are 
particularly  well  represented.  Model  effects  noted  in 
previous  studies  (a  spurious  anticyclone  and  front,  and  a 
displaced  northward  branching  of  the  Gulf  Stream)  distort 
the  mean  flow  in  other  locations.  The  spatial  distribution 
of  mean  kinetic  energy  (MKE)  is  likewise  distorted,  although 
the  frequency  distribution  of  MKE  and  the  total  MKE  over  the 
ocean  appear  to  be  basically  correct.  Mean  flow  vector 
diagrams  show  the  model  fails  to  capture  much  of  the 
disorganized  nature  of  surface  flow  in  the  North  Atlantic. 
This  observation  is  reflected  in  plots  of  standard  deviation 
ellipses,  which  show  the  model's  variance  to  be  low  over 
most  of  the  North  Atlantic,  particularly  so  in  the  south  and 
southeast  Subtropical  Gyre. 

Lagrangian  integral  time  scales  for  the  1/6-degree 
model  are  overestimated,  on  average,  by  factors  of  about 
about  1.9  and  1.8  in  the  zonal  and  meridional  directions, 
respectively.  Integral  length  scales  are  slightly 
underestimated  on  average,  by  factors  of  slightly  greater 
than  unity  and  about  1.2  in  the  zonal  and  meridional 
directions,  respectively. 

The  higher  vertical  and  horizontal  resolution  of  the 

1/10-degree  model  make  a  significant  improvement  in  the 

spatial  distribution  of  energy  and  representation  of  the 

mean  flow.  The  magnitudes  of  variance  are  greatly  improved, 
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the  rich  eddy  structure  of  the  Gulf  Stream  evident,  and  the 
complex  recirculation  patterns  of  the  mid-Subtropical  Gyre 
much  closer  to  that  revealed  by  the  drifters.  The  weak 
Azores  Front  is  resolved  in  the  1/10 -degree  model,  whereas 
it  was  absent  in  the  1/6-degree  model.  The  spurious  features 
of  the  1/6-degree  model  are  either  absent  or  greatly 
attenuated  in  the  1/10-degree  version. 

Our  evaluation  indicates  that  the  displaced  branching 
of  the  North  Atlantic  Current,  spurious  circulation  effects, 
under-resolved  variability,  and  limited  feature  definition 
of  the  1/6-degree  model  make  it  inadequate  for  use  in  a 
global  coupled  forecasting  model.  Although  the  1/10-degree 
comparison  is  considered  preliminary  on  the  basis  of  its 
abbreviated  time  domain,  the  results  of  the  higher 
resolution  run  clearly  indicate  substantial  improvement  in 
the  above  areas  and  are  strong  enough  to  nominate  the 
1/10-degree  model  for  consideration  in  the  coupled  system. 
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I .      INTRODUCTION 

The  paradigm  for  Navy  forecasting  of  ocean  conditions 
in  the  next  millennium  is  that  of  a  global  high-resolution 
atmosphere/ocean/ice  coupled  model  with  the  ability  to 
assimilate  data.  Currently  computer  resources  are 
unavailable  for  such  as  system;  however,  as  multi-processor 
super  computing  technology  continues  to  evolve,  it  is 
appropriate  to  evaluate  candidate  high-resolution  global 
ocean  models  for  potential  use  in  a  coupled  forecasting 
system.  This  is  done  by  quantitatively  assessing  the  ability 
of  the  models  to  reproduce  features  and  processes  with  short 
time  scales  (up  to  2  weeks) .  In  this  thesis  we  evaluate  one 
particular  high-resolution  model,  the  Los  Alamos  National 
Laboratory  (LANL)  Parallel  Ocean  Program  (POP)  model.  Two 
different  formulations  are  considered,  a  1/6-degree  near- 
global  version,  and  a  recent  advancement  of  a  1/10 -degree 
version  covering  only  the  North  Atlantic  Ocean. 

Surface  drifting  buoys  provide  one  of  several  standards 

•by  which  to  assess   the   fidelity  of  general  circulation 

models  in  representing  ocean  flow.  The  frame  of  comparison 

they  offer  is  unique,  as  they  suffer  neither  the  immobility 

of  moored  instruments  nor  the  coarser  temporal  resolution 

(-10   day  resampling  period)   associated  with  altimetry. 

Drifter  data  have  been  archived  for  well  over  two  decades 

and   some   regions,   the   North   Atlantic   Ocean   basin   in 

particular,   have   been   extensively   sampled.   An   improved 

method  for  conditioning  raw  drifter  data,   in  place  since 
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1996  (Hansen  and  Poulain,  1996),  has  produced  a  large, 
consistent  data  set  of  velocity  values.  These  are  provided 
as  uniform  six-hour  time  series,  facilitating  both  Eulerian 
and  Lagrangian  computations.  This  study  compares  output  from 
the  1/6-degree  and  1/10-degree  models  to  data  from  North 
Atlantic  surface  drifters  on  the  basis  of  Eulerian  and 
Lagrangian  velocity  statistics.  • 

Previous  comparisons  of  the  1/6-degree  POP  model  to 
actual  data  sources  include  TOPEX/POSEIDON  (T/P)  altimetry 
(McClean  et  al .  ,  1997),  Pacific  basin  surface  drifters 
(Lemon,  1997),  current  meters  and  altimetry  (T/P  and  GEOSAT) 
(Maltrud  et  al .  ,  1998),  and  tide  gauges  (Tokmakian,  1996). 
The  1/10-degree  model  has  been  compared  to  satellite 
altimetry  data  and  various  in  situ  sources,  and 
intercompared  with  the  1/6-degree  model  (Smith  et  al .  , 
1999)  . 

The  purpose  of  this  study  is  to  provide  an  independent 
assessment  of  the  model's  reproduction  of  the  North  Atlantic 
surface  circulation  in  both  the  1/6-degree  and  1/10-degree 
versions.  The  results  of  this  analysis  will  address  the 
tradeoff  between  adequate  grid  resolution  and  the  need  to 
minimize  machine  time  and  memory;  the  latter  concern  is 
significant  as  we  consider  a  global  coupled  forecasting 
system.  Specific  objectives  of  the  study  are  to  characterize 
regions  of  agreement  and  disagreement  between  model  and 
data,   assign   likely   causes   for   areas   of   significant 


departure,  and  to  determine  locations  where  the  1/6-degree 
formulation  is  adequate. 

Several  intercomparison  schemes  are  used,  based  on  the 
following  data  and  model  output: 


1.  Actual  drifter  velocities  at  each  recorded  track 
point,  referred  to  simply  as  "drifter  data". 

2.  Velocities  obtained  from  numerically  simulated 
(appendix  A)  trajectories  calculated  from  the  1/6- 
degree  model's  velocity  fields.  The  modeled  particles 
are  "released"  at  the  same  locations  and  times  as  their 
real  counterparts.  This  set  of  trajectories  will  be 
referred  to  by  the  name  "same  start". 

3.  Modeled  drifter  velocities,  as  in  (2),  but  using 
random  start  times  and  locations.  There  are  615  random 
deployments,  corresponding  in  quantity  to  the  actual 
drifter  releases.  The  record  lengths  of  these  simulated 
drifters  are  also  randomized,  following  an  exponential 
distribution  with  mean  lifetime  equal  to  that  of  the 
actual  drifters.  This  trajectory  set  will  be  referred 
to  by  the  term  "random  start". 

4.  Velocities  at  model  grid  points  closest  to  actual 
drifter  measurements  in  time  and  space.  The  extraction 
of  these  velocities  was  done  by  applying  a  binary  mask 
to  the  model  lattice,  hence  this  output  will  be 
referred  to  as  "masked" .  Masked  output  sets  are 
generated  using  both  the  1/6-degree  model  grid  and  the 
1/10 -degree  model  grid. 


A  comparison  between  1  and  4  (both  1/6  and  1/10  model 
versions)  will  be  made  on  the  basis  of  Eulerian,  or  location 
based,  statistics.  Lagrangian,  or  trajectory-based, 
statistics  will  be  used  to  compare  1  with  2.  Because  only 
two  years  of  1/10-degree  model  output  were  available  at  the 
time  of  these  analyses,  only  Eulerian  comparisons  are  made 


for  this  version.  A  full  description  of  the  Eulerian  and 
Lagrangian  statistics  is  given  in  Chapter  IV. 

Finally,  a  comparison  between  sets  2  and  3  will  give  a 
rough  estimate  of  the  degree  of  bias  incurred  in  the  1-2 
comparison  through  differences  in  deterministic  and  random 
sampling  patterns. 

The  model /data  comparisons  are  made  quantitatively,  and 
simple  statistical  rejection  tests  for  equal  means  are 
conducted  under  the  null  hypothesis  that  the  model  output 
does  not  significantly  differ  from  the  actual  data.  This 
approach  is  conservative  in  that  it  seeks  to  identify  areas 
of  significant  deficiency  in  the  model's  representation  of 
the  surface  circulation. 

Chapter  II  provides  an  overview  of  the  surface 
circulation  of  the  North  Atlantic  Ocean,  including 
relatively  recent  findings.  It  is  provided  as  a  framework 
for  understanding  our  results  and  those  of  others  cited. 
Chapter  III  gives  a  history  and  description  of  both 
1/6-degree  and  1/10-degree  models,  and  describes  the  drifter 
data  set.  A  review  of  previous  POP  comparisons  is  also 
provided.  Chapter  IV  presents  the  methods  used  in  the 
Eulerian  analysis,  and  describes  our  results.  Chapter  V 
describes  the  methods  and  results  of  the  Lagrangian 
analysis.  Finally,  Chapter  VI  summarizes  our  findings  and 
offers  recommendations. 


II.   SYNOPSIS  OF  CIRCULATION  AND  ITS  VARIABILITY  IN  THE 

NORTH  ATLANTIC  OCEAN 

This  section  is  intended  to  familiarize  the  reader  with 
the  mean  circulation  patterns  and  variability  in  the  North 
Atlantic  ocean,  as  an  aid  to  understanding  the  findings  of 
this  work  and  others  cited.  Figure  1  illustrates  the  major 
currents  referenced  herein.  Figure  2  shows  the  topography  of 
the  Atlantic  Ocean. 

The  basin-wide  mean  circulation  is  driven  by  two  gyres, 
the  North  Atlantic  Subpolar  and  Subtropical  Gyres, 
themselves  driven  by  the  large-scale  wind  stress  curl.  The 
Subpolar  Gyre  consists  of  the  East  Greenland  Current  (EGC) 
and  West  Greenland  Current  (WCG)  to  the  north,  the  Labrador 
Current  to  the  west,  the  North  Atlantic  Current  (NAC)  to  the 
south,  and  the  Irminger  Current  (IC)  to  the  east.  The 
Subtropical  Gyre  is  formed  by  the  Gulf  Stream  to  the  west 
and  northwest,  the  NAC  to  the  northeast,  the  Portugal  and 
Canary  Currents  to  the  west,  the  North  Equatorial  Current 
(NEC)  to  the  south,  and  the  Antilles  Current  and  Florida 
Current  to  the  southwest.  There  are  additionally 
recirculation  currents  within  the  gyre,  which  transport 
water  toward  the  NEC . 

The  two  gyres  meet  in  a  broad  confluence  south  of  the 
Newfoundland  Grand  Banks.  Here  cold,  relatively  fresh  water 
from  the  Labrador  shelf  comes  into  contact  with  the  warmer, 
more  saline  waters  of  the  Sargasso  Sea.  It  is  a  site  of 
intense  eddy  activity,  where  Gulf  Stream  meanders  form  both 
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Figure  1.  Mean  Currents  and  Fronts  of  the  North  Atlantic 
Ocean.  Abbreviations  used  are  for  the  East  Iceland  Current 
(EIC) ,  East  Greenland  Current  (EGC) ,  West  Greenland  Current 
(WGC)  ,  Irminger  Current  (IC),  Antilles  Current  (AC),  Loop 
Current  (LC) ,  Caribbean  Count ercurrent  (CCC) ,  Jan  Mayen 
Front  ( JMF) ,  Norwegian  Current  Front  (NCF) ,  Iceland-Faroe 
Front  (IFF),  Subarctic  Front  (SAF),  and  Azores  Front  (AF)  . 
After  Tomczak  and  Godfrey,  (1994)  . 
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Figure  2.  Bottom  topography  of  the  North  Atlantic  Ocean. 
Bathymetry  is  from  the  ET0P05  database,  courtesy  of  the 
National  Oceanographic  Data  Center. 
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cyclonic  and  anticyclonic  rings,  transporting  cold  shelf 
water  southward  and  warm  Sargasso  water  northward, 
respectively. 

The  greatest  current  velocities  in  the  Atlantic  Ocean 
are  observed  in  the  Gulf  Stream,  specifically  between  Cape 
Hatteras  and  Cape  Cod.  Johns  et  al . (1995)  analyze  the  output 
from  an  array  of  13  current  meters  in  this  area,  finding  a 
maximum  velocity  of  about  2  00  cm/s  at  the  surface, 
decreasing  to  70  cm/s  at  1000  m.  It  is  noteworthy  that  Johns 
et  al .  (1995)  employ  a  "stream-coordinate"  system,  where  the 
coordinate  axes  are  always  aligned  with  the  instantaneous 
direction  of  the  Gulf  Stream.  They  point  out  that  a  simple 
Eulerian  averaging  with  bins  aligned  to  the  standard  axes 
will  give  an  inaccurate  picture  of  the  mean  flow  if  any 
recirculations  or  counterflows  are  present,  whereas  the 
stream-coordinate  method  will  separate  these  effects  from 
the  larger  flow.  Richardson  (1983)  makes  a  related  point 
regarding  averages,  namely,  that  resolving  long-term  mean 
velocities  in  the  Atlantic  requires  enormous  amounts  of 
data.  This  is  because  most  of  the  energy  is  contained  in  the 
eddy  field  rather  than  the  mean  field,  and  because  the 
dominant  time  scales  are  in  the  range  3  0-100  days. 
Richardson's  (1983)  study  of  eddy  kinetic  energy  (EKE)  in 
the  Atlantic  gives  a  peak  value  in  the  Gulf  Stream  of  3000 
cm2/s2,  based  on  a  2°x2°  binning  scheme,  where  large 
amplitude  meanders  begin  to  form  at  37°N,  67°W.  He  further 
identifies  three  high  energy  tongues  associated  with  the 


Stream's  downstream  course,   the  strongest  coinciding  with 
its  flow  around  the  Grand  Banks  into  the  Newfoundland  Basin. 

The  Gulf  Stream  exhibits  pronounced  temporal  and 
spatial  variability.  The  Florida  Current,  which  feeds  into 
the  Gulf  Stream,  shows  significant  seasonal  variability  in 
transport  levels,  with  March  volume  transports  11  Sverdrup 
(Sv)  higher  than  in  November  (Tomczak  and  Godfrey,  1994). 
Larson  and  Bushnell  (1995)  also  demonstrate  decade-scale 
variations  in  Florida  Current  transport.  Spatially,  mean 
transports  increase  from  about  3  0  Sv  off  Miami,  Florida 
(Niiler  and  Richardson,  1973)  to  a  maximum  of  150  Sv  near 
65°W  (Worthington,  1976)  due  to  contributions  from  gyre 
recirculation  and  cold-core  ring  reabsorption.  Thereafter 
the  Gulf  Stream  begins  to  lose  flow  to  Sargasso  Sea 
recirculation  and  ring  formation,  decreasing  to  120  Sv  by 
-the  Tail  of  the  Grand  Banks  (Worthington,  197  6)  .  Finally, 
there  is  substantial  variability  in  the  Gulf  Stream's 
location,  on  both  annual  and  interannual  scales.  In  a  five- 
year  survey,  Auer  (1987)  shows  the  mean  position  between 
70°W  and  44 °W  to  reach  a  northern  extreme  in  September  and  a 
southerly  limit  in  February,  with  the  magnitude  of  the 
annual  shifts  increasing  to  a  maximum  near  the  Newfoundland 
Grand  Banks . 

The  part  of  the  subtropical  gyre  just  east  of  50°W  is 

known  as  the  Gulf  Stream  Extension.  Here  the  North  Atlantic 

Current  (NAC)  breaks  off  as  a  northward  flow,  makes  an  acute 

turn  to  the  southeast,  then  just  as  sharply  veers  to  the 
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northeast.  This  sinuous  course  forms  a  loop  in  the  vicinity 
of  50°N,  42°W  which  has  become  known  as  the  "Northwest 
Corner".  Beyond  this  point,  the  NAC  continues  as  a  frontal 
jet  reaching  velocities  of  100  cm/s  (Krauss,  1986). 
Southward  of  the  jet  a  strong  eddy  field  extends  westward  to 
the  Mid-Atlantic  Ridge.  EKEs  in  this  field,  based  on  a  3°x3° 
grid,  range  from  1500  crnVs2  south  of  Newfoundland  to  600 
cm2/s2  in  the  Northwest  Corner  to  3  00  cm2/s2  at  the  Mid- 
Atlantic  Ridge  (Krauss  and  Kase,  1984). 

As  the  NAC  continues  to  the  northeast  it  loses  flow  to 
the  northward- flowing  Irminger  Current,  which  branches  off 
Southwest  of  Iceland.  It  then  sheds  a  small  amount  of  flow 
to  the  south,  which  supplies  the  Portugal  and  Canary 
currents.  By  this  point  the  NAC  has  slowed  and  diffused  into  . 
a  slow  northeast  drift  in  the  direction  of  Scotland  (Pickard 
and  Emery,  1990)  . 

The  region  to  the  west  and  southwest  of  Iceland  is 
characterized  by  low,  seasonally  variable  velocities 
generally  below  5  cm/s  and  punctuated  by  a  band  of  EKE  of 
about  200  cm2/s2  (based  on  2°  longitude,  1°  latitude  bins) 
coinciding  with  the  Sub-Arctic  Front  (SAF)  (Otto  and  Van 
Aken,  1996) .  The  currents  west  of  Iceland  and  around  the  tip 
of  Greenland  are  primarily  driven  by  thermohaline  gradients, 
with  speeds  of  up  to  7  0  cm/s  observed  in  the  outer  edges  of 
both  the  East  and  West  Greenland  Currents  (Krauss,  1995) . 

The  Norwegian  Atlantic  Current  (NwAtC) ,  also  called  the 

Norwegian  Current,  is  formed  as  the  convergence  of  the  NAC 
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and   the   East   Iceland   Current   (EIC)  ,   with   additional 

contribution  from  a  weak  northeast  drift  to  the  west  of  the 

Iceland-Faroe   Front   (IFF).   Within   a   rectangle   bounded 

roughly  by  Iceland,  Norway,  Spitzbergen,  and  the  east  coast 

of  Greenland  are  the  Greenland,  Iceland,  and  Norwegian  Seas, 

collectively  called  the  Nordic  Seas.  Beneath  these  lie  the 

Norwegian,  Greenland  and  Lofoten  Basins,  and  the  Iceland 

Plateau.   The  seas  are  bounded  on  the  east  by  the  West 

Spitzbergen  Current,  a  branch  of  the  NwAtC,  and  to  the  north 

by  the  EGC.  Poulain  et  al .  (1996),  study  the  tracks  of  107 

drifters  released  in  the  vicinity  of  Iceland  and  the  IFF. 

The  tracks  reflect  earlier  evidence  that  the  NwAtC  passes 

along   the   coast   of   Norway   as   separate   coastal   and 

continental  margin  flows.  They  measure  maximum  velocities  of 

over  110  cm/s  where  the  flow  branches  reassimilate  off  the 

north  end  of  the  Norwegian  peninsula.  The  currents  along  the 

entire  periphery  of  the  seas  appear  to  be  driven  by  sharp 

thermohaline   fronts,   though   the   authors   also   reveal 

topography  to  be  a  major  influence.  In  particular,  Poulain 

et  al .   (1996)   identify  a  western  branch  of   the  NwAtC 

deflected  by  sloping  bathymetry  into  the  region  of  Jan  Mayen 

Island,  and  the  segregation  of  the  larger  circulation  into 

separate  sub-basin  gyres  by  submarine  ridges.  Finally,  their 

study  shows  EKE  to  scale  roughly  with  the  mean  flow  in  this 

region,  with  maximum  values  close  to  700  cm2/s2  (based  on  2° 

longitude,  1°  latitude  bins) ,  corresponding  to  the  maximum 

flow  found  off  the  northwest  Norwegian  coast.  The  Iceland 
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Plateau,   in  contrast,   is  distinctively  low  in  both  mean 
kinetic  energy  (MKE)  and  EKE. 

The  drifter  study  of  Poulain  et  al .  (1996)  represents 
one  of  many  that  have  greatly  aided  our  knowledge  of  the 
North  Atlantic  circulation.  Richardson  (1983),  using  110 
buoys,  maps  EKE  across  the  entire  North  Atlantic  basin, 
revealing  the  rich  turbulent  structure  of  the  Gulf  Stream 
and  noting  that  EKE  measurements  from  ship  drift 
underestimate  the  true  values  by  a  factor  of  about  two.  In  a 
similar  study  of  the  eastern  North  Atlantic,  Krauss  and  Kase 
(1984)  identify  a  homogeneous  pool  of  low  EKE  water  (less 
than  10  0  cm2/s2,  based  on  a  3°x3°  grid)  southward  of  about 
40°N,  in  agreement  with  Richardson's  (1983)  2°x2° 
measurements  of  100  cm2/s2  in  the  eastern  North  Atlantic  and 
North  Equatorial  Current,  and  200  cm2/s2  in  mid-gyre. 

The  works  of  Krauss  (1995)  and  Otto  and  Van  Aken  (1995) 
have  already  been  cited;  additionally,  both  studies  clearly 
reveal  an  northeastward  flow  in  the  region  south  of  Iceland, 
which  divides  at  the  Faroe  Islands  to  enter  the  Norwegian 
Sea  via  the  Iceland-Faroe  Ridge  and  the  Faroe-Shetland 
Channel.  Krauss'  (1995)  drifters  also  show  the  Irminger 
Current  to  split  into  westward  and  eastward  flowing  branches 
at  the  Reykjanes  Ridge.  Both  these  studies  add  valuable 
insight,  as  knowledge  about  the  weak  currents  in  this  region 
is  particularly  poor. 

Colin  de  Verdiere  (1983)   investigates  the  Lagrangian 

properties  of  16  drifters  deployed  in  a  100  km  diameter  area 
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in  the  eastern  North  Atlantic,  finding  good  agreement  up  to 
3  0  days  with  homogeneous  dispersion  theory.  Based  on  centers 
of  gravity  of  energy  spectra,  Lagrangian  time  scales  are 
about  12  days,  and  Eulerian  time  scales  about  2  6  days. 
Krauss  and  Boning  (1987)  study  the  tracks  of  113  drifters  in 
eddy  field  associated  with  the  Gulf  Stream  Extension,  again 
finding  good  agreement  with  homogeneous  dispersion  theory 
over  60  days,  and  a  decreasing  trend  in  Lagrangian  time 
scales  from  the  subtropics  to  the  NAC .  They  find  eddy  length 
scales  to  be  fairly  homogeneous  across  the  Atlantic,  and 
eddy  diffusivity  (the  integral  of  the  Lagrangian 
autocovariance)  to  increase  by  a  factor  of  about  four  from 
30°N  to  50°N. 


13 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


14 


III.  MODEL  OUTPUT  AND  DRIFTER  DATA  SETS 

A.    MODEL  DESCRIPTION,  OUTPUT,  AND  PREVIOUS  COMPARISONS 

The  POP  model  is  a  3 -dimensional ,  z-level,  eddy- 
resolving,  primitive  equation  ocean  general  circulation 
model.  The  code  derives  from  that  developed  by  Bryan  (1969), 
Cox  (1970,  1984)  and  Semtner  (1974) .  Semtner  and  Chervin 
(1988)  modify  the  code  for  a  1/2°  grid  and  improve  the 
realism  of  the  topography.  Further  refinements  produce  a 
1/4°  version  (Semtner,  1997)  which  includes  a  free  surface 
formulation  by  Killworth  et  al .  (1991)  to  allow  the  use  of 
unsmoothed  topography.  These  improvements  increase  the 
model's  EKE  and  topographically-related  eddy  energy  in  the 
Southern  Ocean,  where  comparisons  to  Geosat  altimetry  had 
indicated  deficits  (Wilken  and  Morrow,  1994) . 

The  1/6-degree,  2  0  level,  POP  model  is  a  further 
refinement  of  this  1/4°  version;  it  is  reformulated  for 
massively  parallel  architecture  (Smith  et  al .  ,  1992)  and 
uses  an  implicit  free-surface  method  for  the  barotropic  mode 
(Dukowicz  and  Smith,  1994).  The  spatial  domain  extends  from 
75°S-65°N  to  77°S-77°N.  Salinity  and  temperature  restoring 
to  Levitus  (1982)  climatology  is  implemented  over  the  upper 
kilometer  of  the  water  column  in  the  zone  poleward  of 
latitude  70°  (Maltrud  et  al . ,  1998). 

A  higher  resolution  1/10°,  40-level  simulation  of  the 

North  Atlantic  only  is  run  by  Smith  et  al .   (1999).  The 

spatial  domain  is  20°S  -  72°N,  98°W  -  17°E,  which  includes 

15 


the  Gulf  of  Mexico  and  the  western  Mediterranean  Sea.  The 
1/10°,  40-level  North  Atlantic  run  evaluated  in  this  study 
is  an  advancement  of  the  above  model;  it  includes  the  mixed 
layer  formulation  known  as  K-Profile  Parameterization  (KPP) 
developed  by  Large  et  al .  (1994)  .  This  later  run  is  forced 
with  daily  wind  stresses  from  the  Navy  Operational  Global 
Atmospheric  Prediction  System  (NOGAPS) ,  while  the  former  run 
is  forced  with  wind  stresses  from  the  European  Centre  for 
Medium  Range  Weather  Forecasts  (ECMWF) .  For  clarity  in  the 
following  discussion,  the  latter  will  be  referred  to  as  the 
NOGAPS  version  and  the  former  the  ECMWF  version. 

1.    Forcing  and  Equilibration:  1/6-Degree  Model 

Semtner's  (1997)  1/4-degree  model  is  initialized  using 
the  output  from  a  35  year  integration  at  1/2° ,  (which  in 
turn  is  initialized  from  Levitus  (1982)  climatology),  and  a 
subsequent  5  year  equilibration  at  1/4°.  Following  a  3-year 
run,  the  end  state  is  interpolated  to  1/6°  resolution  and 
used  to  initialize  the  1/6°  version,  which  is  allowed  to 
equilibrate  for  5  years.  (McClean  et  al . ,  1997) 

Three  consecutive  runs  (POP5,  P0P7 ,  POP11)  are  made  at 

1/6°,  with  the  last  three  initialized  from  snapshots  of  the 

end  state  of  their  predecessors.  The  three  sequential  runs 

are  made  to  investigate  the  effects  of  different  wind  and 

heat   forcings.   Salinity   and   temperature   restoration   is 

obtained   by   interpolation   from   Levitus   (1982)   monthly 

climatology,  and  wind  forcing  fields  are  interpolated  from 
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2.5°  wind  grids  provided  by  ECMWF .   POPS  uses  3 -day  wind 
fields  temporally  interpolated  from  ECMWF  monthly  averages, 
while  P0P7  and  P0P11  use  3 -day  fields  averaged  from  twice- 
daily  ECMWF  winds.   P0P11  additionally  employs  heat  flux 
forcing  from  a  method  developed  by  Barnier  et  al .   (1995), 
based  on  ECMWF  climatological  fields.  (Maltrud  et  al .  ,  1998) 
A  fourth  run  is  initialized  from  P0P11  output  for  the 
end  of  1992  and  run  through  1997.  Wind  forcing  is  supplied 
by  daily  ECMWF  wind  stresses  interpolated  to  every  model 
time  step  (30  minutes),  rather  than  by  the  monthly  or  3- 
daily  averages  of  the  previous  runs.  Upper  level  velocities 
for  the  run  under  study  are  saved  as  daily  snapshots. 

2.    1/6-Degree  Model  Grid  and  Bathymetry 

The  1/6-degree  POP  model  is  defined  over  a  lattice  of 

2  0   concentric   Mercator   grids,    one   for   each   depth 

representation.  On  a  given  depth  surface,  there  are  1280 

longitudinally-spaced   grid   points    for   each   of   896 

represented  latitudes.  The  longitude  range  is  global;  the 

latitudes   are   defined   from  77°S   to   77°N.   The  vertical 

spacing  is  25  m  for  the  first  four  levels,   and  varies 

between  3  5  m  and  55  0  m  for  the  remainder.   The  average 

horizontal  resolution  is  -1/5°  (22.42  km  at  the  surface) 

with  lowest  resolution  of  0.28°   (surface,  -31  km)   at  the 

equator   to   0.06°   (surface,   -6.5   km)   at   the   poleward 

boundaries.   Because   the   Rossby   radius   of   deformation 

decreases  poleward  from  the  equator,   a  Mercator  grid  is 
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advantageous  as  its  horizontal  resolution  also  increases 
toward  the  poles.  The  POP  grid  allows  resolution  of  the 
first  baroclinic  Rossby  radius  at  all  latitudes  \X\    <    ~    40°; 

hence,  the  model  should  be  expected  to  represent  eddy 
formation  processes  at  least  up  to  the  mid-latitudes. 
(Maltrud  et  al.  ,  1998) 

The  bathymetry  is  interpolated  to  the  1/6°  grid  from 
the  National  Geophysical  Data  Center's  high-resolution 
(1/12°)  ETOP05  database. 

3.    Forcing  and  Equilibration:  1/10-Degree  Model 

The  ECMWF  1/10-degree  model  is  initialized  from  June 
Levitus  (1982)  temperature  and  salinity  climatology,  and 
allowed  to  spin  up  from  rest,  beginning  1  June  1985,  for  a 
period  of  5.3  years.  Wind  and  salinity  forcing  are  as  per 
the  1/6-degree  version;  the  Barnier  et  al .  (1995)  heat 
fluxes  are  also  used,  except  that  solar  radiation  is  treated 
as  an  independent  quantity  and  allowed  to  penetrate  to 
depth.  At  the  end  of  the  spinup  (1  October  1990),  the  wind 
stress  is  reset  to  1  October  1995  and  the  integration 
continued  until  1  July  1996.  (Smith  et  al.,  1999). 

The  NOGAPS  1/10-degree  run  is  initialized  using  the 
restart  file  from  the  ECMWF  run  at  the  end  of  1992.  This 
version  uses  daily  NOGAPS  wind  stresses  temporally 
interpolated  to  every  model  time  step  (7  minutes)  and 
spatially  interpolated  from  a  1.25°  grid  to  the  model  grid. 


Output  from  the  first  two  years  of  this  1993-1997  run  forms 
the  1/10-degree  database  for  this  study. 

4.  1/10-Degree  Model  Grid  and  Bathymetry 

The  horizontal  sections  of  the  1/10-degree  model 
lattice  are  mercator  grids,  with  resolutions  varying  from 
11.1  km  at  the  equator  to  3.2  km  at  the  northern  boundary. 
In  contrast  to  the  1/6-degree  model,  the  1/10-degree  model 
resolves  the  first  baroclinic  Rossby  radii  at  all  latitudes 
up  to  60 °N.  There  are  40  vertical  levels,  varying  in 
thickness  from  10  m  at  the  surface  to  250  m  at  depth  (Smith 
et  al.,  1999).  Again,  ETOP05  bathymetry  is  interpolated  to 
the  model  grid. 

5.  Output  Fields  and  Binning 

Zonal  (u)  and  meridional  (v)  velocity  fields  from  the 

midpoint  of  the  1/6-degree  model's  upper  level  (12.5  m)  are 

saved  as  daily  snapshots  spanning  the  period  January  1,  1993 

to  December  31,  1997.  Likewise,  daily  second  level  (15.0  m) 

snapshots  from  the  1/10-degree  model  are  saved  for  the 

period  January  1,  1993  to  December  31,  1994.  These  depths 

are  picked  to  most  closely  match  the  depth  of  the  drifters 

(15m) .  The  final  spatial  domain  of  the  analysis  is  chosen  as 

80°W  -  20°E,  20°N  -  70°N,  an  area  defined  by  the  bulk  of  the 

drifter  data.  Both  drifter  data  and  model  output  are  binned 

onto  2°x2°  and  5°x5°  grids  for  the  Eulerian  and  Lagrangian 

comparisons,   respectively,   with   the   exception   that   no 
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Lagrangian  statistics  are  calculated  for  the  1/10-degree 
model . 

6.    Results  of  Previous  Comparative  Studies 

This  section  gives  a  brief  description  of  the  results 
of  previous  comparisons  of  the  POP  model  with  various  data 
sources  in  our  region  of  interest. 

Maltrud  et  al .  (1998)  provide  a  comparative  evaluation 

of  P0P5,  P0P7,  and  P0P11,  and  remark  that  on  the  whole  each 

faithfully   reproduces   the  wind-driven   circulation.   They 

note,   however,   that   problems   remain   in   the   model's 

resolution   of   the   turbulent   energy   spectrum   and   its 

representation  of  meridional  overturning.  Though  the  stream 

functions  from  all  three  simulations  closely  resemble  the 

profiles   of   the   Sverdrup   balance,   where   most   of   the 

departures  can  be  explained  by  thermohaline,  topographic,  or 

inertial   effects,   Maltrud   et   al .    (1998)   point   to   a 

significant  exception  in  the  excessively  poleward  separation 

of  the  northern  hemisphere  western  boundary  currents.   In 

particular,   the  Gulf  Stream  separation  point  at  37°N  is 

displaced  northward  in  the  model  by  a  large,   persistent 

anticyclone,  an  artifact  of  the  model.  The  authors  note  that 

this  effect  may  be  a  native  feature  of  Brian-Cox  type  models 

run  at  medium  to  high  resolutions,  though  recent  runs  at 

both  1/6-degree  using  greater  vertical  resolution  (Chao  et 

al.,   1996)  and  at  1/10-degree   (Smith  et  al .  ,   1999)   show 

improved  realism.  The  addition  of  ECMWF  surface  heat  flux 
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forcing  to  the  formulation  for  temperature  restoration  give 
POP-11  an  improved  mass  transport  in  the  lower  Gulf  Stream, 
though  all  three  versions  of  the  model  suffer  from  an  under- 
represented  proportion  of  flow  through  the  Florida  Strait 

(Maltrud  et  al .  ,  1998).  This  is  apparently  due  to  a  loss, 
versus  augmentation,  of  flow  through  the  Windward  Passage, 
caused   possibly   by   exaggerated   resistance   through   the 
Yucatan  Channel  (Maltrud  et  al .  ,  1998).  Both  the  Atlantic 
and  global  meridional  heat  transports,  and  the  phase  and 
amplitude   of   sea   surface   temperature   and   height,   are 
improved  by  the  use  of  heat  flux  forcing  (Maltrud  et  al .  , 
1998);  difficulties  remain,  however,  in  representing  western 
boundary  undercurrent  depth  and  water  properties,  partly  due 
to  biases   in  initial   conditions.   The  runs  using  3   day 
averaged  wind   fields   (POP-7   and   POP-11)   produce   great 
improvement  in  sea  surface  height  (SSH)  variability,  though 
the  observed  values  are  still  low  compared  to  altimetry 

(Maltrud  et  al . ,  1998). 

McClean  et  al .  (1997)  compare  the  mesoscale  variability 
of  the  1/6-degree  model  output  (along  with  a  1/4-degree 
version  called  the  Parallel  Ocean  Climate  Model  (POCM) )  to 
TOPEX/ POSEIDON  (T/P)  altimetry.  Their  objectives  are  to 
assess  the  impact  of  decreased  friction  and  increased 
resolution  between  the  1/4-degree  and  1/6-degree  models,  and 
to  validate  the  models  on  the  basis  of  SSH  variability,  EKE, 
and  length  scales. 
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The  POP  model  is  found  to  explain  about  60%  of  the 
variability  globally,  and  the  POCM  about  50%.  Locally  the 
POP  outperforms  the  POCM  and  closely  tracks  with  T/P  in  the 
vicinity  of  strong  currents;  otherwise  both  models  generally 
underestimate  variability  by  about  50%.  The  spatial 
distribution  of  the  variability  resembles  that  of  T/P, 
except  for  an  excessive  eastward  displacement  of  the 
Northwest  Corner  (Fu  and  Smith,  1996)  .  In  all  versions  of 
the  model,  EKE  is  found  to  be  underestimated,  due  possibly 
to  the  quality  of  the  wind  forcing  fields,  the  periodicity 
of  the  wind  averaging,  and  inadequate  modeling  of  the  mixed 
layer.  Length  scales  are  more  realistic  in  the  POP  than  in 
the  POCM,  though  both  models  correctly  reproduce  their 
spatial   distributions,  as  observed  by  T/P  altimetry. 

To  infer  the  relative  effects  of  increased  resolution 
over  decreased  friction  between  the  two  models,  McClean  et 
al .  (1997)  calculate  a  zonally-averaged  meridional  section 
of  SSH  in  a  latitudinal  region  where  both  models  resolve  the 
first  baroclinic  Rossby  radius,  and  argue  that  the  observed 
absence  of  a  significant  increase  in  SSH  variability 
indicates  the  friction  coefficient  plays  a  subordinate  role 
(refer  to  Figure  3)  .  This  result  is  supported  by  Treguier 
(1992),  using  a  comparison  of  EKEs  in  a  similar  model,  by 
varying  only  the  mixing  coefficients. 

Tokmakian  (1996)  compares  output  from  the  POCM  and  POP 

to  global  tide  gauge  data,  finding  that  the  quality  of  the 

wind  forcing  fields  is  the  primary  factor  affecting  the 
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Figure    3 .    Meridional    section   of    zonally-averaged   sea   surface 
height    variance     (cm2)     from    POP-11     (solid    line)     and    P0CM-4B 
(dashed   line)    models.    After  McClean   et   al .     (1997). 
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model's  representation  of  SSH.  Specifically,  she  notes  no 
significant  differences  in  the  models'  performance  at 
latitudes  below  35°N.  Fu  and  Smith  (1996)  note  that  sea 
level  variance  is  generally  underestimated  by  a  factor  of 
about  two,  mainly  in  areas  of  high  eddy  activity,  and 
suggest  that  better  spatial  resolution  is  required  to  fully 
resolve  mesoscale  eddies.  They,  too,  find  an  improvement 
when  heat  flux  forcing  is  substituted  for  temperature 
restoration.  Finally,  Fu  and  Smith  (1996)  note  that  the 
model  represents  fairly  well  the  intraseasonal  fluctuations 
at  and  above  midlatitude  and  the  interannual  variability  of 
the  tropics. 

Smith  et  al .  (1999)  compare  results  of  the  1/10-degree 
model  to  various  in  situ  sources  (e.g.,  hydrography, 
drifters,  current  meters),  satellite  data,  and  to  output 
.  from  the  POP-5/POP-7/POP-11  sequence.  They  note  the  1/10- 
degree  model  produces  a  substantial  improvement  over  the 
1/6-degree  model's  representation  of  the  time-averaged  North 
Atlantic  circulation,  and  find  the  former  much  more  closely 
captures  the  shape  and  magnitude  of  the  mesoscale  energy 
spectrum.  Smith  et  al .  (1999)  find  the  1/10-degree  model  to 
distribute  about  70%  of  the  total  basin  kinetic  energy  as 
EKE,  and  note  the  spatial  distribution  of  this  EKE  to  be 
much  more  in  line  with  observations  than  that  of  the  1/6- 
degree  model . 

Smith  et  al .   (1999)  note  continuing  difficulties  in 

modeling  the  Gulf  Stream  and  the  North  Atlantic  Current  at 
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1/10°,  however.  They  find  the  separation  point  of  the  North 
Atlantic  Current,  in  the  region  of  the  Northwest  Corner,  to 
lie  too  far  to  the  northwest.  The  Gulf  Stream  also  appears 
too  far  south,  and  its  meander  envelope  too  wide  to  the  west 
of  the  New  England  Seamounts . 

B.    DRIFTER  DATA  SET 

In  1988  the  World  Climate  Research  Program  organized 
the  Global  Drifter  Program  (GDP)  to  collect,  interpret,  and 
disseminate  sea  surface  temperature  and  velocity  data  from 
globally  deployed  drifters.  One  of  GDP's  objectives  was  to 
design  and  build  a  relatively  inexpensive,  long-lived 
satellite-tracked  drifter  with  specific  water  following  and 
slip  characteristics.  Their  drifter  was  developed  at  the 
Scripps  Institution  of  Oceanography  in  La  Jolla,  California. 
It  is  hereafter  referred  to  as  the  WOCE/TOGA  Lagrangian 
drifter,  as  the  specifications  are  those  of  the  Surface 
Velocity  Program  sponsored  by  those  international  programs. 
Sybrandy  and  Niiler  (1991)  provide  details  of  its 
construction  and  a  parts  list  in  the  WOCE/TOGA  Lagrangian 
Drifter  Construction  Manual.  (Niiler,  1995). 

In  1988  a  coordinated  effort  to  seed  the  tropical 
Pacific  Ocean  with  GDP  drifters  was  begun  as  the  "Pan- 
Pacific  Current  Study".  This  project  was  followed  in  1991 
with  Atlantic  deployments,  which  provide  the  drifter 
database  for  this  study. 
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1.  The  WOCE/TOGA  Lagrangian  Drifter 

The  GDP's  Lagrangian  drifter  consists  of  a 
cylindrically-f ramed  nylon  "holey-sock"  drogue  tethered  to  a 
fiberglass  surface  float.  The  length  of  the  tether  is  such 
that  the  center  of  the  drogue  travels  in  the  mixed  layer  at 
15  m.  Figure  4  shows  the  dimensions  of  the  entire  drifter 
assembly.  Within  the  surface  buoy  are  housed  submergence  and 
sea  surface  temperature  sensors,  and  a  battery,  transmitter 
and  antenna  for  transmitting  temperature  readings  to 
satellite-carried  service  Argos  instruments.  Due  to  the  cost 
of  satellite  connectivity  (Gianetti,  1993)  the  transmitters 
are  programmed  to  be  on  for  one  day,  then  off  for  two  days. 
The  battery  packs  are  designed  for  two  years  life  under  this 
duty  cycle.  The  submergence  sensor  is  installed  as  an  aid  to 
detecting  drogue  loss,  as  the  floats  frequently  submerge 
when  drogued  (Niiler,  1995)  . 

2 .  Data  Collection  and  Conditioning 

Drifter  transmitter  signals  are  received  passively  by 

ARGOS  instruments  aboard  polar-orbiting  National  Oceanic  and 

Atmospheric  Administration  (NOAA)  satellites.  There  are  at 

least  two   such  satellites   operating,   in  sun- synchronous 

orbits,  at  a  given  time.  Their  equatorial  crossings  occur  in 

a  fixed  westward  procession  at  102 -minute  intervals.  Mean 

pass  duration  is  ten  minutes,  and  the  repetition  period  for 

a  drifter's  transmissions  is  90  seconds.   (Service  Argos, 

1984) 
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Figure  4 .  Diagram  of  the  WOCE-TOGA  Global  Lagrangian 
Drifter,  showing  dimensions  in  cm.  After  Poulain  et  al . 
(1996)  . 
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Each  transmission  localizes  the  drifter,  via  doppler 
shift,  to  a  cone  formed  with  the  satellite  at  the  apex.  The 
orbital  elements  of  the  satellite,  based  on  its  reception  of 
reference  signals  from  fixed  ground  transmitters,  further 
restrict  the  possible  locations  to  the  branch  of  a 
hyperbola.  The  intersection  of  the  hyperbolas  formed  from 
the  satellite  pass  finally  place  the  drifter  at  one  of  two 
possible  locations.  Argos  then  resolves  the  final  ambiguity 
using  geof easibility  criteria. 

Position  and  temperature  data  are  then  forwarded  to 

NOAA's  Atlantic  Oceanographic  and  Meteorological  Laboratory 

(AOML)  .  The  data  received  from  Argos  are  subject  to  errors 

from   instrument   instability,   variations   in   number   of 

transmissions   received,   and   inaccuracies   in   measuring 

tracking  geometry;  the  temporal  spacing  is  uneven  as  well 

(Hansen  and  Poulain,   1996)  .   To  make  the  raw  data  more 

suitable   for   analysis,   AOML   applies   an   editing   and 

interpolation  technique  developed  by  Hansen  and  Poulain 

(1996),  and  briefly  outlined  below,  to  remove  any  erroneous 

values  and  create  uniform  six-hour  time  series  of  velocity 

and  temperature. 

Bad  position  data  are  identified  using  a  sequential 

point-rejection  algorithm,  where  a  new  position  is  rejected 

if  the  velocity  it  implies  compared  to  the  previous  position 

exceeds  a  certain  threshold.  The  threshold  is  specified  from 

known  circulation  characteristics  of  the  region  or  from 

histograms   of   measured   velocity.   Temperature   data   are 
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considered  bad  if  they  exceed  a  running  weighted  mean,  where 
the  weights  have  been  determined  experimentally.  For  both 
positional  and  temperature  data,  the  respective  algorithms 
are  applied  in  the  forward  and  reverse  directions.  Only 
values  that  are  offenders  in  both  directions  are  rejected, 
so  that  good  data  are  not  unnecessarily  lost.  (Hansen  and 
Poulain,  1996) 

Once  the  position  and  temperature  data  are  screened  for 
errors,  they  are  optimally  interpolated  to  uniform  6-hour 
time  series.  Position  data  are  treated  as  separate  latitude 
and  longitude  time  series;  velocity  values  are  then  obtained 
by  center-differencing  the  interpolated  longitude  and 
latitude  values  (Poulain,  1993). 

The  interpolation  method  is  a  unidimensional  variant  of 
kriging,  a  procedure  used  in  mining  and  hydrology,  where 
interpolated  values  are  calculated  as  linear  combinations  of 
neighboring  measurements .  The  weights  are  normally 
determined  by  minimizing  the  kriging  variance,  or  mean 
squared  difference  between  true  and  interpolated  values. 
Because  the  interpolation  variable  in  the  present 
application  is  position,  however,  the  kriging  variance  as 
defined  is  not  guaranteed  to  have  general  existence.  It  is 
therefore  specified  in  terms  of  an  analytic  structure 
function  model  whose  parameters  are  determined  empirically 
using  actual  drifter  data.  (Hansen  and  Poulain,  1996) 
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3 .    Drogue  Slippage 

A  major  concern  in  the  design  of  the  WOCE/TOGA 
Lagrangian  drifter  was  that  it  have  a  predictable  slip 
through  the  water  (Sybrandy  and  Niiler,  1991) ,  where  slip  is 
defined  as  the  difference  between  drogue  velocity  and  the 
current  velocity  averaged  over  the  length  of  the  drogue. 
Additionally,  to  avoid  compounding  error,  it  was  desired  to 
keep  the  slip  within  the  limits  of  the  average  daily 
velocity  uncertainty  incurred  from  ARGOS  tracking.  This 
uncertainty  was  estimated  at  1  cm/s,  based  on  a  conservative 
value  of  1  km  for  positional  error  (Niiler  et  al . ,  1995)  . 

Nonlinear  motion  effects  were  controlled  by  choosing  a 
spherical  float  shape  to  minimize  aliasing  of  wave-induced 
motion  into  the  horizontal  plane,  a  thin,  stiff  tether  to 
reduce  looping  or  kicking  from  waves,  and  a  drogue  shape 
that  was  resistant  to  kiting.  To  investigate  the  water- 
following  characteristics  of  the  resulting  design,  a  series 
of  field  tests  was  performed  under  various  water  shear  and 
wind  conditions  using  a  drifter  outfitted  with  calibrated 
current  meters,  one  at  each  end  of  the  drogue.  The  average 
meter  velocity  was  taken  as  the  slip  speed  (Sybrandy  and 
Niiler,  1991) . 

From  these  tests,  wind  is  deduced  to  be  the  dominant 
force  producing  slip,  and  a  statistical  model  is  offered  to 
correct  for  velocity  measurements  using  the  formula 
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U=3.5(W/R)  +  9.5(DU/R).  (1) 

where  U  (cm/s)  represents  the  slip  velocity,  W  (m/sec)  the 
wind  velocity  at  10  m,  DU  (cm/s)  the  velocity  difference 
across  the  drogue,  and  where  R  is  defined  as  the  ratio  of 
the  product  of  the  drogue  area  and  its  drag  coefficient  to 
the  sum  of  the  products  of  the  tether  and  submerged  floats 
and  their  drag  coefficients.  The  equation  associates  a  slip 
speed  of  1  cm/s  and  a  wind  speed  of  10  m/sec  with  an  R  value 
of  3  5;  the  WOCE/TOGA  drifter  has  an  R  value  of  40,  so  that 
the  design  criterion  is  met.  (Niiler,  1995) 

4.    Wind  Correction 

Niiler  et  al .  (1995)  provide  a  detailed  description  of 
the  above  experiment,  but  add  that  no  buoy  deployments  were 
made  in  winds  above  10  cm/s  because  the  floats  proved 
undeployable  in  the  resulting  sea  conditions.  Subsequent 
regression  analysis  of  a  much  more  complete  worldwide 
database  of  drifters  and  global  wind  products  has  produced  a 
simpler  but  improved  formula  for  the  WOCE/TOGA  drifter 
(Stephen  Pazan,  personal  communication) : 

U       =U  -0.12W  (2) 


corr     uncorr 


where  U        is  in  cm/s,  U  is  in  cm/s,  and  W   is  in  m/s.  We 

corr  '      uncorr  ' 

use  the  latter  equation  in  this  study  to  compute  corrected 
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velocities  for  both  Eulerian  and  Lagrangian  analyses.  The 
multiplier  W  is  obtained  for  each  observation  by  spatial 
bilinear  interpolation  of  the  corresponding  daily  NOGAPS 
wind  components . 

5.    Temporal  and  Spatial  Domains 

Our  drifter  data  set  includes  624  of  the  GDP  drifters 
deployed  in  the  Atlantic.  The  total  record  length  extends 
from  January  1,  1993  to  December  31,  1997.  The  drifter  data 
are  subsampled  at  daily  intervals  (every  4  data  points)  to 
match  the  daily  sampling  of  model  output.  Figure  5  shows  the 
temporal  distribution  of  the  subsampled  drifter 
observations.  The  number  of  deployed  drifters  is  between  50 
and  100  for  nearly  every  day  in  the  time  domain.  Significant 
peaks  occur  in  October,  1993  (127  drifters),  and  March  of 
1997  (126  drifters)  .  Years  1994  and  1995  see  nearly 
monotonic  decreases  and  increases,  respectively,  with  the 
lowest  population  in  the  time  domain  (46  drifters)  occurring 
in  January,  1995.  Drifter  population  is  near  100  for  most  of 
1996. 

The  spatial  domain  is  chosen  to  include  the  greater 
part  of  the  drifter  coverage.  Though  the  Gulf  of  Mexico, 
western  Caribbean,  and  latitudes  as  far  south  as  12°S  are 
represented  by  the  drifters  in  this  5  year  period,  the 
coverage  is  considered  too  sparse  for  meaningful 
comparisons.  Additionally,  we  choose  70°N  as  an  upper 
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Figure  5.  Temporal  distribution  of  daily  drifter 
observations,  January  1,  1993  to  December  31,  1997,  over  the 
geographical  region  20°N  -  70°N,  80°W  -  20°E. 
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boundary  to  remain  outside  the  restoring  buffer  zone.  Our 
choice  of  an  eastern  limit  is  driven  simply  by  the  northern 
boundary  and  land  constraints.  We  define  our  region  of 
interest  as  80°W  to  20°E,  20°N  to  70°N. 

Figure  6  gives  the  deployment  locations  of  the  actual 
drifters,  while  Figure  7  shows  the  "deployment"  locations  of 
the  simulated  drifters  for  the  "random  start"  scheme.  There 
is  no  discernable  pattern,  naturally,  for  the  "random 
starts".  The  actual  drifter  plot,  however,  is  punctuated  by 
concentrated  releases  off  Cape  Cod  and  the  coasts  of  Iceland 
and  Norway,  between  the  Azores  and  Canary  Islands,  east  of 
the  Mid  Atlantic  Ridge  between  44 °N  and  49 °N,  and  in  the 
Labrador  Sea.  In  mid-ocean,  the  deployments  are  otherwise 
fairly  uniformly  scattered. 

6.    Drifter  Tracks  and  Simulated  Trajectories 

The  "spaghetti"  plot  of  actual  drifter  tracks,  Figure 
8,  reflects  the  dispersion  from  the  deployment  sites,  with 
high  track  densities  corresponding  to  the  areas  of  intense 
deployment.  The  coherent  Gulf  Stream  and  northern  boundary 
currents  are  clearly  seen,  while  the  southeast  basin  is 
characterized  by  circuitous,  rambling  patterns.  An  eastward 
surface  flow  across  the  Iceland-Faroe  Ridge  is  clearly  shown 
by  the  concentrated  deployments  south  of  Iceland,  whereas 
the  path  of  the  North  Atlantic  Current  through  the  Rockall 
Channel  into  the  Faroe-Shetland  Channel  is  sparsely  sampled. 
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Figure  6.  Deployment  sites  of  actual  drifters  released  in 
geographical  region  20°N  -  70°N,  80°W  -  20°E,  during  time 
period  January  1,  1993  to  December  31,  1997. 
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Figure  7.  Sites  of  615  simulated  drifter  "deployments"  from 
the  "same  start"  regime.  Spatial  and  temporal  boundaries  are 
those  of  the  actual  drifters,   20°N  -  70°N,   80°W  -  20°E, 
January  1,  1993  to  December  31,  1997. 
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Figure  8.  "Spaghetti"  diagram  of  actual  drifter  tracks, 
January  1,  19  93  to  December  31,  19  97.  Sampling  frequency  is 
daily. 

37 


Notable  "holes"  in  the  coverage  occur  over  the  eastern 
boundary  of  the  Labrador  Sea,  the  Iceland  Plateau  and  region 
north  of  Iceland,  the  wide  continental  shelves  off  Canada 
and  the  southern  North  American  coast,  in  the  marginal  seas 
off  western  Europe,  and  southwest  of  the  Canary  Islands. 

Tracks  from  the  "random  start"  regime  (Figure  9) 
reflect  more  frequent  inner  shelf  "deployments",  especially 
in  the  Labrador  Sea.  There,  the  ring  of  coherent  boundary 
flow  is  in  stark  contrast  to  the  more  static  central 
Labrador  Sea.  The  simulated  drifters  released  west  of 
Scotland  correctly  show  the  parallel  northeast  flows  of  the 
North  Atlantic  Current  in  the  Rockall  Channel  and  the  drift 
currents  north  of  the  Rockall  Plateau.  The  flow  across  the 
Iceland-Faroe  Ridge  agrees  with  that  indicated  by  the 
drifters . 

Coverage  is  lower  in  other  areas,  such  as  the  Iceland 

shelf,  which  is  nearly  devoid  of  drifters  except  to  the 

north,   and  the  region  west  of  the  Strait  of  Gibraltar. 

Significantly,   the   absence   of   "random   start"   tracks 

outlining   the   Reykjanes   Ridge   reflects   Krauss '   (1995) 

bifurcation   of   the   Irminger   Current.   Other   lapses   in 

coverage   occur   in   the   same  places   as   for   the   actual 

drifters.  The  northwest  Nordic  Seas,  southern  Newfoundland 

Grand  Banks,  and  the  eastern  reaches  of  the  Labrador  Sea  are 

sparsely  covered.  Except  for  the  North  Sea,   the  European 

marginal  seas  are  infrequently  sampled  by  the  "random  start" 

trajectories . 
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Figure  9.  "Spaghetti"  diagram  of  numerical  drifter  tracks 
from  the  "random  start"  simulation.  January  1,  1993  to 
December  31,  1997.  Sampling  frequency  is  daily. 
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Track  density  in  the  Subtropical  Gyre  from  roughly  3  0°N 
to  40°N  is  lower  for  the  model  than  for  the  actual  drifters. 
The  more  direct  "random  start"  trajectories  are  due  to  the 
model's  overly  uniform  gyre  recirculation  patterns,  as 
described  in  section  IV. B. 2.  A  comparison  of  the  drifter 
tracks  with  the  "same  starts"  trajectories  (Figure  10)  more 
clearly  reflects  the  differences  in  flow  dynamics.  South  of 
45°N  the  model  has  two  very  strong  zonal  signals,  one 
associated  with  the  Gulf  Stream,  and  the  other  with  a 
spurious  front  to  be  described  in  section  IV. B. 2.  In 
contrast,  the  actual  drifter  tracks  have  a  more  homogeneous 
distribution.  Another  notable  difference  concerns  the 
northward  turn  of  the  Gulf  Stream  along  the  Newfoundland 
Grand  Banks.  Whereas  the  data  show  the  flow  continuing 
northward  into  the  Northwest  Corner,  the  "same  start"  tracks 
show  a  premature  eastward  banking  at  roughly  45  °N, 
consistent  with  the  previously  observed  eastward 
displacement  of  this  feature. 

Figure  11  shows  the  spatial  distribution  of  the  drifter 

data  in  number  of  daily  observations  per  2°x2°  bin.  The 

peaks  in  the  regions  of  concentrated  deployments  noted  above 

are  obvious.  The  entrainment  of  the  drifters  into  the  Gulf 

Stream,   Irminger   Current,   East   Iceland   Current,   and 

Norwegian  Atlantic  Current  is  evident,   as  is  the  broad 

westward  diffusion  of  drifters  released  in  the  southeast 

Subtropical  Gyre.  The  stasis  in  the  central  Labrador  Sea  is 

apparent  as  well. 
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Figure  10.  "Spaghetti"  diagram  of  numerical  drifter  tracks 
from  the  "same  start"  simulation.  January  1,  1993  to 
December  31,  1997.  Sampling  frequency  is  daily. 
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Figure  11.  Number  of  daily  observations  of  actual  drifters 
per  2°x2°  bin,  January  1,  1993  to  December  31,  1997. 
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IV.   EULERIAN  ANALYSIS 

A.    FORMULAE  AND  METHODS 

For  the  Eulerian  analyses  we  partition  our  geographical 
domain  (section  III.B.5)  into  50  zonal  bins  and  25 
meridional  bins.  For  each  bin  we  extract  velocities  from  the 
midnight  drifter  observations  taken  within  its  boundaries, 
and  the  model  velocities  at  the  corresponding  time-location 
coordinates.  We  then  compute  the  Eulerian  statistics 
described  below  for  both  sets  of  velocity  values.  All 
observations  from  the  actual  drifters  are  corrected  for 
drogue  slippage  as  described  in  section  III.B.4. 

Separate  comparisons  are  made  for  the  1/6-degree  model 
run  and  the  1/10-degree  run.  Because  five  years  of  output 
were  available  from  the  former,  whereas  only  two  years  of 
output  were  available  from  latter,  we  use  two  separate 
drifter  data  sets  for  consistency  in  the  comparisons.  The 
drifter  data  for  the  1/6-degree  comparison  is  the  full  data 
set  spanning  the  period  January  1,  1993  to  December  31, 
1997,  while  that  for  the  1/10-degree  comparison  is  an 
abbreviated  set  running  from  January  1,  1993  to  December  31, 
1994. 

For  the  1/6-degree  model  comparisons,   the  data  and 

model  velocities  are  extracted  over  the  time  domain  January 

1,   1993  to  December  31,   1997.  For  the  1/10-degree  model 

comparisons  the  velocities  are  collected  for  the  period 

January  1,  1993  to  December  31,  1997. 
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If  zonal  velocity  (cm/s)  is  denoted  u,     and  meridional 
velocity  (cm/s)  as  v,  then  we  may  define  mean  kinetic  energy 
(MKE,    cia/s2)  ,    by 


MKE^2{(uyE+(V)l)  (3) 


where  the  operator  <>E  denotes  a  simple  average  taken  over 
all  the  data  falling  into  a  given  bin,  and  where  the  prime 
indicates  the  residual  about  the  Eulerian  mean.  Eddy  kinetic 
energy  {EKE,    cmVs2)  is  defined  by 


EKE  =  -{\+^2)  (4) 


where  Xx  and  Z2  (\  >  Z2)  are  the  principal  components  of 
variance,  computed  (Freeland  et  al . ,  1975)  as  the  roots  of 

^)s-X)^)e-4-{u'V%=0.  (5, 

The  roots  A1     and  A2     are  therefore  the  eigenvalues  of  the 

Eulerian  sample  covariance  matrix.  The  two  eigenvalues  and 
their  corresponding  eigenvectors  define  an  ellipse 
describing  the  distribution  of  velocity  variance  for  the 
given  bin;  equivalently,  a  standard  deviation  ellipse  is 
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defined.  The  rotation  angle  6   for  the  ellipse  is  calculated 
(Freeland  et  al . ,    1975)  via 


( 1  ^         (uVO 
-  an(2fl)  =  ,2\       f         .  (6) 

v2J  iv  )E-{U  h 


The  choice  of  bin  size  is  made  to  ensure  adequate 
feature  resolution  and  number  of  observations  per  bin.  It 
should  be  noted  that  the  chosen  resolution  has  an  effect  on 
the  measured  proportion  of  EKE  and  MKE.  Figure  12,  a  plot 
of  EKE  and  MKE  vs.  resolution,  illustrates  this  effect.  The 
points  on  the  curves  are  generated  by  averaging  the  drifter 
EKE  and  MKE  values  for  all  bins,  for  each  given  bin  size. 
Binning  schemes  are  optimal  with  respect  to  robustness  at 
resolutions  where  the  two  curves  run  nearly  horizontally; 
results  are  then  equally  valid  for  any  resolution  in  their 
plateaus.  This  is  unfortunately  not  the  case  for  high 
resolutions,  due  to  the  increasingly  steep  ramping  of  MKE  as 
bin  size  decreases  below  about  six  degrees.  Our  chosen  bin 
size  (2°x2°)  lies  between  these  two  extremes,  enabling  good 
feature  resolution  while  avoiding  the  region  where 
sensitivity  increases  quite  sharply. 

Plots  of  mean  vectors,  principal  standard  deviation 
ellipses,  and  eddy  and  mean  kinetic  energies  for  both  model 
and  drifter  data  are  provided  in  the  subsequent  section. 
Additionally,   we  test  for  significant  departures  of  the 
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model  from  the  "true"  mean  drifter  velocities.  This  is  done, 
bin  by  bin,  using  both  a  "graphical"  method,  and  Hotelling's 
T2  statistic  for  bivariate  normal  populations. 


EKE  and  MKE  vs.  Bin  Size 
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Figure  12.  Variation  of  eddy  kinetic  energy  (EKE)  (crnVs2) 
and  mean  kinetic  energy  (MKE)  (cm2/s2)  with  bin  size 
(degrees) .  The  points  on  the  curves  are  obtained  by  binning 
the  actual  daily  drifter  data  into  bins  of  the  given  size. 
The  present  study  uses  2°x2°  bins  for  the  Eulerian  analysis. 
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For  both  methods  we  extract  subsets  of  the  model  and 
data  velocities  that  are  considered  temporally  and  spatially 
independent .  These  we  denote  by 


Y.  ■[%.*-, V.J  ■  and  vd=  ev  v V.3' 


where  the  subscripts  m  and  d  refer  to  the  model  and  data, 

respectively,   nm  and  nd  are   the   number   of   independent 

observations  from  the  model  and  data,   respectively,   and 

where  the  boldface  indicates  the  vector  quantity  of  both 

zonal   and   meridional   components .   To   obtain   the   above 

subsets,  we  first  subsample  the  daily  velocity  records  at  a 

time   interval   greater   than   a   representative   average 

Lagrangian  integral  time  scale  for  the  North  Atlantic  Ocean, 

taken  as  two  days  (see  section  V.B.I)  .  For  the  observations 

corresponding  to  each  independent  time  step,   a  heuristic 

network   algorithm   is   applied   to   select   a   spatially 

independent   sample.   The   algorithm   first   identifies   all 

observations  that  are  within  a  representative  mesoscale  eddy 

radius  for  the  North  Atlantic.  This  is  taken  to  be  100  km, 

based  on  a  midlatitude  average  of  the  first  internal  Rossby 

radius  of  deformation.  It  then  randomly  chooses  a  point  to 

reject,   reevaluates   the  remaining  network  for  distances 

between  points,  and  repeats  until  no  point  is  within  the 

threshold  distance  of  any  of  its  neighbors.  The  resulting 

points  form  a  candidate  set;  if  the  algorithm  is  applied 

47 


multiple  times,  the  largest  set  among  the  candidates  is 
taken  as  the  sample  for  that  time  step.  The  union  of  all 
such  sets  over  all  time  steps  then  forms  a  temporally  and 
spatially  independent  sample  for  that  bin.  We  apply  the 
network  algorithm  10  times  for  each  independent  time  step, 
over  each  bin. 

In  the  graphical  test  for  equality  of  means,  we  first 
form  95%  standard  error  ellipses  for  the  model  and  data  by 
scaling  the  principal  component  eigenvalues  of  Vm  and  Vd, 
respectively,  by 


t  a      J—  and   t  a 


where  ta  n  is  the  1-a  quantile  of  Student's  t  distribution 

with  n  degrees  of  freedom.  For  each  bin  we  connect  the  tails 
of  the  model  velocity  vector  and  the  drifter  velocity 
vector.  The  standard  error  ellipses  are  then  placed  with 
their  centers  at  the  tips  of  their  respective  vectors.  If 
the  ellipses  fail  to  intersect,  the  velocities  are  assumed 
to  be  significantly  different  for  that  bin. 

Hotelling's  T2  test  (Seber,  1984)  involves  calculating, 
for  V  and  V,,  the  test  statistic 
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Ti  :    _^^_(?ffl  -v,)'S;'(vm  -  v„),         where 


nm+nd 


i    "»  i    nrf 

(nm+nd-2) 

and  where  S   and  S,  are  the  model  and  data  covariance 

m  a 

matrices,   respectively,   for  the  bin.  The  null  hypothesis 
that  the  model  velocity  does  not  depart  significantly  from 

the  drifter  velocity  is  rejected  at  level  a  (.05  for  95% 
confidence)  if 


,     2{nm+n.—2)    „  ~ 

rp2  >  m  g  pg  _     <p2 

71     +71     —J  m      d  m      d 

jn  a 


The  test  is  fairly  insensitive  to  departures  from 
normality  and  unequal  covariance  matrices  when  the  two 
sample  sizes  are  equal  (Seber,  1984)  .  We  do  not  meet  this 
condition  in  general,  due  to  edge  effects  in  binning  model 
velocities  (the  data  domain  is  continuous,  while  the  model 
domain  is  discrete) ,  but  the  sizes  are  close  enough  that  the 
test  should  work  reasonably  well. 

Finally,  we  look  for  the  dependence  of  the  model's 
performance  on  energy  level  (both  EKE  and  MKE)  by  plotting 
the  ratios  of  model  kinetic  energy  to  drifter  kinetic  energy 
against  drifter  kinetic  energy   (see   following  section) . 
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Plots  for  the  Eulerian  results  are  provided  in  the  following 
section;  plot  output  is  presented  only  for  bins  with  at 
least  40  observations  for  the  1/6-degree  comparisons,  and  10 
observations  for  the  1/10-degree  comparisons. 

B.    EULERIAN  RESULTS:  1/6-DEGREE  MODEL 

1.    Mean  Flow  and  Eddy  Kinetic  Energy  Indicated  by  the 
Drifters 

Figure  13  shows  the  mean  circulation  derived  from  the 
drifter  data,  over  the  2°x2°  grid,  for  1993  -  1997.  The 
major  circulation  patterns  discussed  in  chapter  I  are 
evident.  The  Gulf  Stream,  East  and  West  Greenland  Currents, 
Labrador  Current,  Irminger  Current,  East  Iceland  Current, 
and  Norwegian  Atlantic  Current  are  particularly  well- 
defined,  as  is  the  North  Atlantic  Current  west  of  25°W.  The 
Northwest  Corner  is  quite  prominent.  The  eastern  limb  of  the 
Subtropical  Gyre  is  also  evident,  the  Canary  Current  and 
North  Equatorial  Current  clearly  visible,  though  not  nearly 
as  strong  as  the  above  currents.  A  division  of  circulation 
in  the  southern  Nordic  Seas  into  two  sub-gyres  is  also 
evident,  with  an  apparent  boundary  along  the  prime  meridian. 

Maximum  drifter  speed  seen  in  the  Gulf  Stream  (22  6 
cm/s)  agrees  well  with  the  2  00  cm/s  measured  by  Johns  et  al 
(1995) .  The  drifters  reach  85  cm/s  in  the  Norwegian  Atlantic 
Current  (119  cm/s  in  the  6-hourly  observations),  roughly  in 
agreement  with  the  110  cm/s  noted  by  Poulain  et  al .  (1996)  . 
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Peak  speeds  observed  in  the  East  and  West  Greenland 
Currents  (123  cm/s,  142  cm/s)  are  notably  higher  than  the  70 
cm/s  reported  by  Krauss  (1995),  due  perhaps  to  a  more 
extensive  data  set  and  sampling  closer  to  the  cores  of  the 
currents.  Mean  values  for  the  region  south  of  Iceland  and 
west  of  Scotland,  most  of  which  are  from  0.5  to  5  cm/s,  are 
on  the  same  order  as  the  0.7  cm/s  to  4.7  cm/s  cited  by  Otto 
and  Van  Aken  (1996) . 

A  plot  of  EKE  for  the  drifter  data  (not  shown)  reveals 

the  three  tongues  of  high  EKE  noted  by  Richardson  (1983); 

these   are   evident   as   well   in   the   principal   standard 

deviation  ellipse  plots  (Figure  18)  .  Eddy  kinetic  energy  in 

the  Gulf  Stream  at  37°N,   67°W  is   2140   cm2/s2  from  the 

drifters,  two-thirds  of  the  value  calculated  by  Richardson 

(1983).   Our  maximum  observed  EKE  of  2957  cm2/s2,   nearly 

identical  with  Richardson's,  occurred  further  to  the  east, 

in   the   bin   defined   by   38°N-40°N,   62°W-64°W.   Values 

downstream  south  of  Newfoundland  agree  with  the  1500  cm2/s2 

found  by  Krauss  and  Kase  (1984);  our  drifters  indicate  over 

twice  the  60  0  cm2/s2  they  cite  for  the  Northwest  Corner, 

however.  East  of  the  Mid  Atlantic  Ridge  our  estimates  (3  00 

cm2/s2)  are  again  in  agreement.  It  should  be  noted  that  the 

most  of  the  drifters  of  Krauss  and  Kase  (1984)  were  drogued 

significantly  deeper  (100  m)  than  ours,  which  may  explain 

the  large  difference  in  Northwest  Corner  estimates.  Along 

the  Subarctic  Front   southwest  of   Iceland,   EKE  from  the 

drifters  is  in  agreement  with  the  100  cm2/s2  to  200  cm  /s 
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found  by  Otto  and  Van  Aken  (1996)  .  We  find  the  same  low 
values  of  EKE  (<  100  cm2/s2)  noted  by  Richardson  (1983)  and 
Krauss  and  Kase  (1984)  in  the  eastern  basin,  only  bisected 
by  the  third  of  Richardson's  "tongues"  extending  eastward  to 
the  African  shelf.  Our  values  of  EKE  in  the  North  Equatorial 
Current  match  Richardson's  100  cm2/s2  only  over  the  east 
basin;  to  the  west  they  are  about  twice  this  value.  Finally, 
for  the  Nordic  Seas  region  south  of  7  0°N,  we  observe  the 
same  low  (<  100  cm2/s2)  energies  over  the  southern  Iceland 
Plateau  and  higher  values  (200  -  400  cm2/s2)  along  the  course 
of  the  Norwegian  Atlantic  Current,  in  agreement  with  Poulain 
et  al.  (1996)  . 

2 .    Mean  Flow  Comparison 

Figure  14  shows  the  mean  circulation  from  the  masked 
output.  There  are  both  pronounced  similarities  and 
differences  when  compared  to  Figure  13  .  The  unrealistic 
large  anticyclone  noted  by  Maltrud  et  al .  (1998),  where  the 
Gulf  Stream  separates  from  the  continental  shelf  at  37°N,  is 
present  in  the  output.  The  eastward  displacement  of  the 
Northwest  Corner  to  mid-basin,  noted  by  Fu  and  Smith  (1996), 
appears  in  Figure  14,  though  the  northward  turn  of  the  Gulf 
Stream  around  the  Newfoundland  Grand  Banks  is  also 
indicated. 

The  model  reproduces  the  northeastern  drift  through  the 

gap   between   Iceland   and   Scotland.   Though   it   slightly 

exaggerates  its  magnitude,  it  captures  the  lack  of 
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uniformity  in  this  flow.  Equally  chaotic  is  the  flow  on  the 
westward  side  of  the  model's  Northwest  Corner;  here  the 
feature  interferes  with  the  eastern  boundary  of  the  Irminger 
Current.  The  northeast  quarter  of  the  Subtropical  Gyre  is 
well  represented  by  the  model,  but  there  is  very  little 
evidence  of  the  Canary  Current  to  the  south.  The  Azores 
Front  is  absent  as  well  in  the  model  mean  vector  plot,  but 
is  evident  between  32  2N  and  34 2N  in  the  drifter  plot. 

The  model  shows  more  gyre  recirculation  from  the  North 
Equatorial  Current  than  do  the  drifters.  This  may  be  a 
sampling  effect,  however,  due  to  relatively  low  drifter 
concentration  south  of  24°N.  In  the  model,  recirculations 
from  both  the  North  Equatorial  Current  and  the  Gulf  Stream 
veer  eastward  to  form  a  basin-wide  flow  between  2  6°N  and 
30°N.  Smith  et  al .  (1999)  note  the  presence  of  this  effect 
in  the  1/6°  simulations  of  Maltrud  et  al  (1998)  .  According 
to  Milliff  et  al.  (1996),  this  flow  is  the  effect  of  a  front 
in  the  model's  Sverdrup  stream  function  created  by  strong 
wind  stress  curls  off  the  coast  of  Africa  around  30°N. 
Neither  this  feature  nor  the  strong  parallel  patterns  of  the 
model's  recirculation  is  revealed  by  the  drifters. 

In  a  general  sense,  there  is  more  uniform  structure  to 

the  current  patterns  in  the  model  than  is  evidenced  by  the 

drifters.  The  complexity  of  flow  in  the  Gulf  Stream  is 

particularly  under-represented;   the  model   shows  a  wide, 

excessively  uniform  eastward  flow  between  40°N  and  52°N, 

extending  to  the  Mid-Atlantic  Ridge.  Over  this  course  the 
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boundary  between  the  Gulf  Stream  and  the  North  Atlantic 
Current  is  ill-defined  in  the  model,  whereas  the  drifter 
data  show  the  eastward  mean  flow  of  the  Gulf  Stream 
dissipating  west  of  approximately  3  8°W. 

Mean  flow  patterns  and  magnitudes  are  well  reproduced 
above  60°N.  There  is  a  large  degree  of  similarity  between 
model  output  and  drifter  data  along  the  Norwegian  Atlantic, 
East  Iceland,  and  East  and  West  Greenland  Currents.  The 
Labrador  Current  is  also  accurately  represented  down  to 
45°N.  The  model  also  captures  the  separation  of  Nordic  Seas 
circulation  into  sub-gyres  gyres  along  the  ridge  between  the 
Norwegian  and  Lofoten  Basins. 

3 .    Mean  Kinetic  Energy  Comparison 

A  comparison  of  model  and  drifter  MKEs  via  Figures  15 
and  16  shows  differences  in  energy  distribution  along  the 
course  of  the  Gulf  Stream  and  effects  connected  with  the 
model's  artifacts.  There  is  excessive  energy  in  the  Antilles 
Current,  and  along  the  mid-Atlantic  shelf  between  34 °N  and 
40°N,  corresponding  to  the  model's  spurious  anticyclone. 
Along  the  mid-course  of  the  Gulf  Stream,  from  its  seaward 
turning  at  Cape  Hatteras  to  the  Northwest  Corner,  MKEs  are 
underestimated.  The  model's  displacement  of  the  Northwest 
Corner  creates  artificially  low  MKEs  in  the  feature's  actual 
location.  In  contrast,  the  model's  delayed  northward  turn  of 
the  Gulf  Stream  and  its  shifted  Northwest  corner  cause  MKE 
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Figure  15.  Mean  kinetic  energy  (MKE)  (cm  /s  )  for  the  actual 
drifter  data  over  the  5-year  time  domain,  January  1,  1993  to 
December  31,  1997.  Bins  with  fewer  than  40  daily 
observations  contain  no  information  and  are  in  white. 
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Figure  16.  Mean  kinetic  energy  (MKE)  (cm /s  )  for  the  masked 
1/6-degree  model  output  over  the  5 -year  time  domain,  January 
1,  1993  to  December  31,  1997.  Bins  with  fewer  than  40  daily 
observations  contain  no  information  and  are  in  white. 
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to  be  exaggerated  in  the  regions  28°W-38°W,  38°N-46°N  and 
26°W-36°W,  54°N-60N°,  respectively. 

The  energy  associated  with  the  North  Atlantic  Current 
is  too  high  west  of  24°W  and  north  of  52°N.  Together,  these 
effects  indicate  that  though  there  is  obvious  spatial 
mismatch  in  MKE,  the  total  over  the  path  of  the  Gulf  Stream 
may  be  reasonable. 

The  ratio  plot  of  model  MKE/drifter  MKE  (not  shown) 
indicates  the  model  exaggerates  MKE  in  the  latitude  band 
26°N  to  30eN,  and  underestimates  MKE  at  latitudes  30°N  to 
34°N.  These  effects  correspond  to  the  model's  artificial 
front  at  the  former  latitudes,  and  the  Azores  Front 
represented  by  the  drifters  at  the  latter.  The  Eastern 
Labrador  Sea  is  a  region  where  MKEs  are  exaggerated  to  a 
large  degree,  in  the  ratio  sense.  The  actual  magnitude  of 
■the  difference  is  quite  small,  however. 

The  MKE  plots  and  MKE  ratio  plots  both  indicate  good 
agreement  between  model  and  data  along  continental  boundary 
currents,  with  the  exceptions  noted  above.  There  is  fairly 
good  agreement  along  the  North  Equatorial  Current  in  the 
eastern  basin. 

Figure  17  displays  histograms  of  drifter  MKE  and  model 
MKE.  The  histograms  reveal  nearly  identical  energy 
distributions,  indicating  that  the  differences  between  model 
and  drifter  MKEs  are  mainly  in  the  spatial  patterns. 

The  bottom  plot  of  Figure  17  shows  the  distribution  of 

model  MKE/drifter  MKE  with  respect  to  drifter  MKE.  A 
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Figure  17.  Plots  (a) , (b) :  Normalized  frequency  histograms  of 
MKE  (cm2/s2)  for  drifter  data  and  masked  1/6-degree  model 
output.  Plot  (c)  :  Variation  of  the  ratio  model  MKE/drif ter 
MKE  with  drifter  MKE  (cm2/s2)  .  Vertical  scale  is  on  a  log 
(base  10)  axis.  For  all  three  plots,  bins  with  fewer  than  40 
daily  observations  for  either  drifter  data  or  model  output 
are  not  represented. 
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significant  vertical  trend  in  the  graph  would  indicate  that 
the  model's  pattern  of  departure  (overestimation, 
underestimation)  from  the  actual  data  had  some  dependence  on 
the  energy  level.  As  indicated  by  the  plot,  there  is  some 
evidence,  though  not  strong,  that  the  model  tends  toward 
underestimation  with  increasing  MKE. 

4.    Comparison  of  Principal  Variance  Components  and 

Eddy  Kinetic  Energy- 
Figures  18  and  19  give  the  principal  standard  deviation 
ellipses  for  the  drifter  and  masked  data.  A  salient  feature 
is  the  model's  underestimation  of  variance  magnitudes  over 
most  of  the  data  domain.  Agreement  is  fairly  close,  however, 
south  and  southeast  of  Greenland,  north  of  Iceland,  east  of 
the  Mid  Atlantic  Ridge  along  the  course  of  the  North 
Atlantic  Current,  and  along  the  Norwegian  Atlantic  Current. 
As  with  MKE,  an  artificial  contrast  in  variational  energy  is 
driven  by  the  displaced  Northwest  Corner  and  delayed 
northward  turn  of  the  Gulf  Stream.  Model  EKEs  are  slightly 
higher  than  drifter  values  in  the  region  26°W-34°W,  42°N- 
56°N;  in  the  actual  Northwest  Corner  they  are 
underestimated.  Model  EKEs  are  exaggerated  in  the  region 
70°W-76°W,  34°N-40°N,  due  at  least  in  part  to  the  spurious 
anticyclone . 

A  plot  of  drifter  EKE  (not  shown)  clearly  shows  the 
Azores  Front  extending  across  the  eastern  basin,  from  4  6°W 
to  the  North  African  coast,  in  a  roughly  4  degree  band.  This 
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feature  is  evident  in  the  ellipse  plot  as  well.  A  similar 
zonal  band  of  relatively  high  EKE  is  indicated  from  24°N, 
apparently  associated  with  the  North  Equatorial  Current's 
contribution  to  the  Sargasso  Sea  recirculation.  The  absence 
of  this  band  in  the  masked  data  is  consistent  with  the 
model's  more  uniform  recirculation  pattern.  Though  not 
evident  in  the  ellipse  plots,  the  model's  spurious  front 
shows  up  in  the  ratio  plot  of  model  EKE/drifter  EKE  (Figure 
20)  .  Energies  are  still  underestimated,  but  the  artificial 
contrast  is  evident  in  the  latitude  band  2  6°N-3  0°N.  With  the 
exceptions  of  the  preceding  features,  the  differences  in 
variational  energy  between  model  and  data  appear  to  be 
mainly  in  scale,  the  spatial  trends  being  otherwise  similar. 
The  histograms  of  Figure  21  reveal  a  unimodal 
distribution  to  the  drifter  EKE,  with  a  most  frequent  value 
of  about  100  cm2/s2.  The  model  EKE  has  a  nearly  monotonic 
distribution,  with  almost  all  values  occurring  below  250 
cm2/s2.  A  histogram  of  the  ratio  model  EKE/drifter  EKE  (not 
shown)  indicates  the  model  under-represents  EKE  in  well  over 
three-quarters  of  the  bins.  Most  values  fall  within  a  range 
of  0.15  to  0.5.  The  bottom  plot  of  Figure  21  shows  that  the 
quantity  model  EKE/drifter  EKE  shows  no  pronounced  trend 
with  respect  to  energy  level. 

5.    Test  for  equality  of  mean  flow 

The   two  tests   for  equality  of  means   described   in 

section  A  are  done  ten  times  for  each  bin  having  at  least  40 
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Figure  20.  Ratio  of  1/6-degree  model  EKE  to  drifter  EKE  for 
the  5-year  time  domain,  January  1,  1993  to  December  31, 
1997.  Bins  with  fewer  than  40  daily  observations  for  either 
drifter  data  or  model  output  contain  no  information  and  are 
in  white. 
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Figure  21.  Plots  (a), (b) :  Normalized  frequency  histograms  of 
EKE  (cmVs2)  for  drifter  data  and  masked  1/6-degree  model 
output.  Plot  (c)  :  Variation  of  the  ratio  model  EKE/drifter 
EKE  with  drifter  MKE  (cm2/s2)  .  Vertical  scale  is  on  a  log 
(base  10)  axis.  For  all  three  plots,  bins  with  fewer  than  40 
daily  observations  for  either  drifter  data  or  model  output 
are  not  represented. 
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observations.   As   described  above,   the   tests   take   into 

account  the  variance  structures  of  the  model  output  and 

drifter  data.  Figure  22  shows  the  results  of  the  tests.  The 

gray  shades  indicate  the  number  of  times  out  of  ten  the 

given  bin  failed  the  given  test  for  equality  of  means.  There 

is  a  high  degree  of  similarity  in  the  results  of  both  tests, 

and  the  salient  characteristics  are  the  same.  Both  indicate 

nearly   binary   results;   generally,    either   significant 

disagreement  is  found  in  almost  all  ten  tests,  or  else  few 

rejections  occurred  at  5%  level  of  significance.  The  spatial 

patterns  are  equally  distinct,   and  reflect  the  earlier 

findings  of  agreement  (more  properly,  lack  of  disagreement) 

in  the  northern  marginal  currents .  The  northern  extremes  of 

the  Gulf  Stream  south  of  Newfoundland,  and  the  northeast  and 

southwest  boundaries  of  the  Subtropical  Gyre  also  see  few 

rejections . 

Significant  rejection  occurs,  however,  where  the  Gulf 

Stream  separates  from  the  continental  shelf,  and  where  flow 

continues  around  the  Newfoundland  Grand  Banks   into  the 

Northwest  Corner.  These  areas  have  already  been  noted  as 

showing   MKE   mismatch;   there   is   significant   direction 

mismatch  west  of  the  Grand  Banks  as  well,  as  a  result  of  the 

model's  delayed  northern  turn.  To  the  south  and  southwest  of 

Iceland  many  rejections  occur  due  to  the  Northwest  Corner's 

artificial  presence.  The  eastern  boundary  of  the  Subtropical 

Gyre  shows  a  large  number  of  rejections.  This  is  consistent 

with  the  lack  of  resolution  of  the  Canary  Current  and  Azores 
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Figure  22.  Results  of  tests  (section  IV. B. 5)  for  equality  of 
mean  velocities  between  drifter  data  and  1/6-degree  model 
output.  Ten  tests  per  bin  were  conducted  to  detect 
significant  departures  from  equality  at  confidence  level  95% 
(a  =  0.05)  .  Each  used  unique  independent  subsets  of  data 
and  output.  The  color  scale  indicates  the  number  of  failed 
tests . 
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Front  in  the  model,  while  the  features  are  clearly  indicated 
by  the  drifters.  The  large  block  of  rejections  over  the 
central  Subtropical  Gyre  south  of  roughly  3  6°N  seems  to  be 
largely  due  to  direction  mismatch.  Here,  the  model's  surface 
transport  has  a  moderate  meridional  component  due  to 
recirculation  from  the  Gulf  Stream  and  North  Equatorial 
Current,  while  the  effect  is  not  as  pronounced  in  the  data. 

C.    EULERIAN  RESULTS:  1/10 -DEGREE  MODEL 

The  following  results  derived  from  1/6-  and  1/10-degree 
model  output  over  the  period  January  1,  1993  to  December  31, 
1994,  and  the  corresponding  drifter  data  spanning  this  time 
domain.  As  the  data  and  model  output  are  fewer  in  number 
than  for  the  1993  -  1997  comparisons,  these  results  should 
be  considered  preliminary. 

1.    Mean  Flow  Comparison 

Figure  23  shows  the  mean  circulation  indicated  by  the 

drifter  data  for  years  1993  -  1994.  Figures  24  and  25  give 

the  corresponding  two-year  plots  for  the  1/6-degree  and 

1/10-degree  models,  respectively.  The  patterns  evident  in 

the  1993  -  1997  drifter  plot  are  still  present  in  the  1993  - 

1994  subset,  though  much  resolution  of  the  Canary  Current  is 

lost.   A  comparison  of   Figure   25   with  Figure   24   shows 

striking  differences.   The  strong,   coherent  recirculation 

patterns   in  the  middle  Subtropical  Gyre  and  the  wide, 

laminar  eastward  flow  between  40 -N   and  52 °N   are  no  longer 
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present.  The  spurious  front  between  2  6-N  and  3  09N  is  still 
indicated,  but  is  much  attenuated  and  of  lesser  latitudinal 
extent  in  the  1/10-degree  model.  There  is  little  direct 
evidence  that  the  1/6-degree  model's  artificial  anticyclone 
has  been  removed,  due  to  the  sparsity  of  data  over  the  shelf 
region  between  Cape  Cod  and  Cape  Hatteras.  However, 
ancillary  analyses  of  mean  sea  level  height  show  no 
anticyclone  at  this  location  (personal  communication,  Julie 
McClean) .  The  Northwest  Corner  now  occurs  in  the  correct 
longitude  band.  It  no  longer  appears  to  interfere  with  the 
Irminger  Current,  and  the  exaggerated  flow  energy  it  created 
to  the  east  of  the  Mid  Atlantic  Ridge  is  gone.  The  weak 
northeast  drift  south  of  Iceland  is  closer  in  magnitude  to 
the  drifters'  values.  The  Azores  Front  now  clearly  appears 
in  the  model  between  32  2N  and  34 SN,  a  result  also  noted  by 
Smith  et  al .  (1999) . 

The  model's  velocity  field  is  too  strong  over  most  of 
the  Gulf  Stream  and  in  the  Northwest  Corner,  however. 
Elevated  velocities  associated  with  the  feature  extend  too 
far  to  the  north  and  west  as  well.  In  contrast,  the  Canary 
Current  still  appears  absent,  though  this  result  must  be 
weighted  with  the  sparsity  of  data  and  the  weakness  of  the 
actual  feature. 

Noting  these  exceptions,  we  see  that  the  1/10-degree 

formulation   eliminates   significant   shortcomings   in   the 

1/6-degree  model.  The  spurious  effects  seen  in  the  latter 

are  either  eliminated  or  greatly  attenuated.  The  1/10-degree 
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model  appears  to  capture  much  more  of  the  turbulent  nature 
of  the  North  Atlantic  Ocean,  while  resolving  most  of  the 
features  indicated  by  the  drifters. 

2 .    Mean  Kinetic  Energy  Comparison 

Plots  of  MKE  (Figures  26  and  27)  reflect  the 
improvements  noted  above;  additionally  MKEs  in  the  North 
Equatorial  Current  appear  to  be  better  matched  by  the  1/10- 
degree  model.  Though  the  model's  amplification  of  MKE  along 
the  Gulf  Stream  is  clear,  the  high  values  associated  with 
the  northward  turn  do  not  extend  as  far  to  the  southeast  as 
the  drifters  indicate.  The  1/10-degree  model  appears  to 
underestimate  MKE  in  the  eastern  reaches  of  the  Gulf  Stream, 
20°W-34°W,  44°N-52°N. 

A  ratio  plot  of  model  MKE/drifter  MKE  (not  shown) 
reveals  the  large  improvement  in  spatial  distribution  of  MKE 
over  the  central  Subtropical  Gyre.  The  strong  zonal  banding 
produced  by  the  model's  spurious  front  and  its  apparent  lack 
of  an  Azores  Front  is  no  longer  present. 

The  histograms  of  Figure  2  8  have  nearly  the  same  shapes 
as  those  for  the  1/6-degree  comparisons,  again  showing 
differences  to  be  mainly  spatial.  The  one  exception  is  that 
MKEs  above  500  cra2/s2  are  roughly  five  times  more  frequent  in 
the  model  than  in  the  drifter  data,  whereas  the  1/6-degree 
model  bins  had  about  two  times  fewer  MKE  measurements  in 
this  band  than  the  drifters.  It  appears  the  1/10-degree 

74 


-3- 

CD 


CO 
CO 


■a 

Q 


E 
o 

UJ 


I 


^ 


4      - ••   i 


T— 


■ 


en      _o 


^^•■.-. 


<  '"'.i- 


■i 


I      \ 


aw  is         s  i 

*'  ni   - 


■ 


,■>,    * 


f  .V:      M    •        ,      f 

■  - 


9     m((  ;■ 


,*-.  S\ 


v.. 

/'  ^ 

■ 

'." 

1 

v 
\ 

■ 

• 

( 

■■•.   r/ 
4  1 

_  5 


CD 
O 
CO 


O 
CO 


O 
CD 


O 
■3- 


O 
C\J 


O 
O 


O 
CO 


O 
CO 


o 


o 

CM 


Figure  26.  Mean  kinetic  energy  (MKE)  (cm /s  )  for  the  actual 
drifter  data  over  the  2-year  time  domain,  January  1,  1993  to 
December  31,  1994.  Bins  with  fewer  than  10  daily 
observations  contain  no  information  and  are  in  white. 
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Figure  27.  Mean  kinetic  energy  (MKE)  (cm /s2)  for  the  masked 
1/10 -degree   model   output   over   the   2 -year   time   domain, 
January  1,  1993  to  December  31,  1994.  Bins  with  fewer  than 
10  daily  observations  contain  no  information  and  are  in 
white . 
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Figure  28.  Normalized  frequency  histograms  of  MKE  (cm2/s2) 
for  drifter  data  and  masked  1/10-degree  model  output.  The 
time  domain  is  January  1,  1993  to  December  31,  1994.  Bins 
with  fewer  than  ten  daily  observations  for  either  drifter 
data  or  model  output  are  not  represented. 
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model  tends  to  slightly  overestimate  MKE  on  the  average, 
whereas  no  strong  pattern  was  shown  by  the  1/6-degree  model. 
This  observation  is  reflected  by  the  histogram  (not  shown) 
of  model  MKE/drifter  MKE,  which  gives  a  shape  similar  to  the 
corresponding  histogram  of  Figure  17,  but  shifted  to  the 
right . 

3.    Comparison  of  Principal  Variance  Components  and 
Eddy  Kinetic  Energy 

Figure  29  shows  the  principal  standard  deviation 
ellipses  for  the  drifter  data  over  the  period  1993  -  1994. 
Figures  3  0  and  31  give  the  corresponding  plots  for  the  1/6- 
degree  model  and  1/10-degree  model,  respectively,  over  that 
time  interval.  The  increase  in  variance  magnitude  achieved 
by  the  1/10-degree  model  is  evidently  quite  significant. 
Resolution  of  variability  over  most  of  the  Subtropical  Gyre 
is  highly  improved.  The  Azores  Front  is  now  clearly  visible, 
where  it  was  not  indicated  in  the  plots  for  the  1/6-degree 
model.  Variability  over  most  of  the  latitudes  north  50°N  is 
well  represented.  The  model  continues  to  underestimate  in 
the  Nordic  seas  north  of  66°N,  though  to  a  smaller  degree. 
The  variability  over  the  Labrador  Basin  is  now  quite  close 
to  that  indicated  by  the  drifters.  Over  the  Rockall  Plateau 
and  west  to  roughly  26°W,  variability  is  underestimated, 
though  only  slightly.  This  information  is  significant 
because  the  1/6-degree  model  variances  for  this  region  were 
masked  by  the  effects  of  the  displaced  Northwest  Corner. 
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The  decrease  in  variance  magnitudes  away  from  the  center  of 
the  Northwest  Corner  is  more  prolonged  in  the  1/10-degree 
model  than  in  the  data.  The  result  is  that  variances  are 
exaggerated  around  its  periphery,  as  evident  in  Figures  29 
and  31.  Variance  is  exaggerated  as  well  in  the  lower  Gulf 
Stream,  where  it  separates  from  the  continental  shelf.  This 
was  the  case  as  well  for  the  1/6-degree  model,  but  the 
overestimation  is  greater  in  magnitude  and  eastward  extent 
in  the  1/10-degree  model.  Elevated  EKEs  (the  mean  of  the 
principal  variance  components)  extend  too  far  south  in  this 
region  as  well,  a  result  noted  by  Smith  et  al .  (1999) . 

Figure  3  2  shows  histograms  of  model  EKE  and  drifter  EKE 
for  the  1/10-degree  model.  A  comparison  with  the  1/6-degree 
model  histogram  of  Figure  21  shows  the  1/10-degree  model 
distributes  much  less  energy  in  the  0-50  cm2/s2  band.  The 
1/10-degree  model  shows  a  notable  amount  of  energy  in  excess 
of  2  000  cm2/s2,  whereas  the  maximum  EKE  seen  in  the  1/6- 
degree  model  was  1641.3  cm2/s2.  The  energy  distribution  of 
the  1/10-degree  model  corresponds  more  closely  to  that  of 
the  drifters  in  the  0  -  1000  cm2/s2  band.  Energies  beyond 
1000  cm2/s2  are  more  prevalent  in  the  1/10-degree  model 
output  than  in  the  drifter  data,  whereas  the  opposite  is 
true  for  the  1/6-degree  model.  In  particular,  the  high 
energies  beyond  2000  cm2/s2  are  not  seen  in  the  drifter  data 
subset.  Apparently  much  of  the  change  in  EKE  between  the 
1/6-degree  and  1/10-degree  formulation  takes  place  at  the 

higher  and  lower  energy  levels.  A  histogram  of  the  ratio 
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Figure  32.  Normalized  frequency  histograms  of  EKE  (cmVs2) 
for  drifter  data  and  masked  1/10-degree  model  output.  The 
time  domain  is  January  1,  1993  to  December  31,  1994.  Bins 
with  fewer  than  ten  daily  observations  for  either  drifter 
data  or  model  output  are  not  represented. 
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model  EKE/drifter  EKE  (not  shown)  reveals  the  1/10-degree 
model  under-represents  EKE  significantly  less  than  does  the 
1/6-degree  model.   The   1/10-degree  histogram's  center  of 
gravity  is  much  closer  to  unity,  whereas  the  center  for  the 
1/6-degree  histogram  was  slightly  less  than  -1/3. 
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V.    LAGRANGIAN  ANALYSIS 

For  the  Lagrangian  analysis  we  divide  the  geographical 
domain  (section  III.B.5)  into  5°x5°  bins,  giving  20  and  10 
bins,  respectively,  in  the  zonal  and  meridional  directions. 
For  each  bin,  we  obtain  drifter  velocities  as  in  the 
Eulerian  analysis,  and  generate  velocities  for  both  "same 
start"  and  "random  start"  regimes  by  the  method  described  in 
Appendix  B.  The  time  domain  for  the  Lagrangian  analysis  is 
the  five  year  period  January  01  1993  to  December  31,  1997. 

A  Lagrangian  description  of  motion  is  based  upon  the 
trajectories  of  individual  particles.  Consequently, 
statistics  based  on  this  approach  are  defined  over  a  time 
series  of  observations.  The  kernel  of  all  such  statistics  is 
the  Lagrangian  autocovariance  function.  Through  Taylor's 
(1921)  theory  of  homogeneous  turbulence,  this  function  can 
be  related  to  the  steady-state  single  particle  eddy 
diffusivity  and  in  turn  to  the  Lagrangian  integral  time  and 
length  scales.  Formal  definitions  of  these  quantities  follow 
below.  The  notation  is  adapted  from  Davis  (1991)  and  Poulain 
et  al.  (1996) . 

A.    FORMULAE  AND  METHODS 

The  Lagrangian  description  requires  a  geo-temporal 
reference,  or  location-time  "tag".  Let  the  Lagrangian 
velocity  and  displacement  at  time  t  of  a  particle  passing 
through  location  x  at  time  t0    be  denoted  by  v(t|x,  t0)  and 
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r(t|x,  t0),  respectively.  Let  the  Lagrangian  operator  <#>L 
signify  an  average  taken  over  an  ensemble  of  particles 
released  randomly  from  x.  Then  we  can  define  the  Lagrangian 
mean  velocity  and  Lagrangian  mean  displacement  at  time 
t0  +  T  of  an  ensemble  of  particles  released  from  x  at  time  t0 
by 

V(r,f0,x)  =  (v(f0+rlx,f0))  (7) 

R(r,r0,x)=(r(f0+Tlx,f0))  (8) 

Finally,  we  denote  the  residual  velocity  and  displacement 
about  these  means  as  v'  and  r' ,  respectively. 

The   Lagrangian   autocovariance   function   is   defined 
(Davis,  1991)  as 

Pl} (x, f)  =  (v,  (r0  I x, t0 )v,  (t0+r\ x, t0 )\  .  (9) 

where  T     is   the   time   lag   of   interest,   and  where   the 

subscripts  i  and  j  may  take  the  values  1  for  the  zonal 
direction  or  2  for  the  meridional  direction.  Davis  (1991) 
further  defines  the  single  particle  eddy  diffusivity  as 

o 
KlJ(x,T)  =  -(v'i(t0\x,t0)r'j(t0-T\x,t0))L  =  jPv(x,T')df.  (10) 
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The  four  components  are  also  frequently  written  in  tensor 
notation  .- 


k(x,t)  = 


jcu(x,t)    ku(x,t) 
_k2](x,t)    tc22(x,r) 


(11) 


Taylor  (1921)  showed  that  for  homogeneous  (location 
independent)  and  stationary  (time  independent)  turbulence 
fields,   the   time   lag   dependence   of  k{x,t)       disappears 

eventually,  as  the  diffusion  becomes  a  random-walk  (Colin  de 
Verdiere,  1983),  leaving 


fC(X,T)  = 


jp} ,  (xr')dr'  Jp]2  (xr')dr' 

-r  -r 

0  0 

JP2]  (xj')dT  JP22  {yL,T)dz' 


K2\        ^22 


=  fC°°,T>T, 


(12) 


where  T  represents  the  time  at  which  the  random  walk  begins. 
At  time  lags  in  excess  of  T,  the  diagonal  elements  of  K  are 
simply  the  products  of  the  respective  velocity  variances  and 
integral  time  scales  (Freeland  et  al . ,  1976).  Equivalently, 


U  = 


(v;(r0lx,r0)2)L|"V,7,  i  =  \,2 
(v;(r0lx,r0)2)^/C  i  =  l,2 


(13) 
(14) 
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give  the  Lagrangian  time  and  length  scales  (Poulain  et  al . , 
1996)  in  terms  of  the  steady-state  diffusivity. 

Our   method   for   calculating   the   sample   Lagrangian 
autocovariance  function  is  presented  in  Appendix  A,  and  is 
applied  to  each  5°x5°  bin.   We  make  the  assumptions  of 
homogeneity  and  stationarity  in  order  to  employ  result  (12) 
above.  Homogeneity  should  clearly  not  be  expected,  however, 
in  regions  where  the  mean  flow  is  significantly  upset  by 
dynamic  or  topographic  effects.   In  our  data,   pronounced 
inhomogeneities   occur   mainly   over   the   western   basin, 
associated  with  the  Gulf  Stream  and  Northwest  Corner.  Lesser 
degrees  of  inhomogeneity  are  generally  evident  along  the 
paths  of  the  other  continental  boundary  currents. 

B.     LAGRANGIAN  RESULTS 

Lagrangian  time  and  length  scales  are  calculated  using 
equations  12  and  13,  respectively,  of  section  V.A.  The 
required  eddy  dif fusivities  are  computed  by  the  method  given 
in  section  V.B.4  below.  Figures  3  3  through  3  6  show  computed 
Lagrangian  time  scales  (zonal  and  meridional)  for  the 
drifter  data  and  "same  starts"  regimes.  Results  are  shown 
for  all  bins  having  at  least  250  observations.  Figures  37 
and  3  8  give  the  corresponding  ratios  of  model  time  scale  to 
drifter  time  scale,  where  results  are  shown  for  bins  having 
at  least  250  observations  for  both  model  and  drifters. 
Figures  3  9  and  40  present  frequency  histograms  of  integral 
time  scales  for  the  drifters  and  both  "same  start"  and 


"random  start"  numerical  trajectories.  Figures  41  through  48 
show  the  above  plots  for  integral  length  scales . 

1.    Integral  Time  Scales 

A  striking  feature  in  all  the  time  (and  length)  scale 
plots  is  the  demarcation  of  the  Gulf  Stream.  In  the  time 
scale  plots,  particularly  those  of  the  meridional  component, 
it  produces  a  swath  of  low  values,  with  longer  scales  to  the 
northwest  and  southeast.  Beyond  this  effect,  time  scales 
tend  generally  to  decrease  with  latitude,  so  that  the 
highest  are  seen  below  the  Gulf  Stream.  The  model 
overestimates  zonal  time  scales  by  factors  up  to  3.6  in  the 
latitude  band  20°N  to  30°N  and  by  up  to  4.9  north  of  the 
Gulf  Stream  and  its  extension.  Outside  these  regions,  zonal 
time  scales  are  in  fairly  good  agreement,  as  shown  by  the 
ratio  plot  of  Figure  37. 

Meridional  time  scales  are  also  exaggerated  north  of 
the  Gulf  Stream,  (by  up  to  5.6  times)  .  There  is  also 
overestimation  in  the  2  0°N  to  3  0°N  band,  but  only  in  the 
eastern  basin,  by  up  to  a  factor  of  3.3.  The  meridional 
ratio  plot  reflects  the  broader  extent  of  values  near  unity, 
particularly  distinct  over  the  North  American  Basin. 

"Random  start"  (not  shown)  and  "same  start"  time  scales 

appear  to  match  fairly  closely  in  magnitude  and  spatial 

distribution.   The  histograms  of  Figures   3  9  and  4  0   show 

similar   frequency  distributions  as  well.   For   the   zonal 

direction,  both  regimes  show  distributions  with  most  time 
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scales  between  four  to  nine  days  and  long,  relatively 
uniform  tails  out  to  15  days.  Neither  shows  time  scales 
below  about  two  days.  In  contrast,  the  drifter  histogram 
indicates  a  significant  number  of  those  low  values,  and  none 
greater  than  eight  days.  For  the  meridional  direction,  most 
model  time  scales  fall  in  the  2-6  day  range  while  higher 
values  occur  out  to  ten  and  twelve  days  for  the  "same  start" 
and  "random  start"  schemes,  respectively.  The  drifter 
histogram  again  shows  a  large  number  of  occurrences  below 
two  days,  none  greater  than  about  eight  days,  and  most 
scales  in  the  one  to  five  day  range.  Clearly,  the  model 
inflates  integral  time  scales  on  the  whole,  with  the  effect 
most  pronounced  in  the  zonal  direction. 
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Figure  33.  Lagrangian  integral  time  scale  (days),  zonal 
component,  drifters.  The  time  domain  is  January  1,  1993  to 
December  31,  1997.  Bins  with  fewer  than  250  daily 
observations  contain  no  information  and  are  in  white. 
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Figure  34.  Lagrangian  integral  time  scale  (days),  zonal 
component,  "same  start"  simulation.  The  time  domain  is 
January  1,  1993  to  December  31,  1997.  Bins  with  fewer  than 
250  daily  observations  contain  no  information  and  are  in 
white . 
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Figure  35.  Lagrangian  integral  time  scale  (days),  meridional 
component,  drifters.  The  time  domain  is  January  1,  1993  to 
December  31,  1997.  Bins  with  fewer  than  250  daily 
observations  contain  no  information  and  are  in  white. 
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Figure  36.  Lagrangian  integral  time  scale  (days),  meridional 
component,  "same  start"  simulation.  The  time  domain  is 
January  1,  1993  to  December  31,  1997.  Bins  with  fewer  than 
250  daily  observations  contain  no  information  and  are  in 
white . 
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Figure  37.  Ratio  of  "same  start"  Lagrangian  integral  time 
scale  to  drifter  integral  time  scale,  zonal  direction.  The 
time  domain  is  January  1,  1993  to  December  31,  1997.  Bins 
with  fewer  than  2  50  daily  observations  for  either  real  or 
simulated  trajectories  contain  no  information  and  are  in 
white. 
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Figure  38.  Ratio  of  "same  start"  Lagrangian  integral  time 
scale  to  drifter  integral  time  scale,  meridional  direction. 
The  time  domain  is  January  1,  1993  to  December  31,  1997. 
Bins  with  fewer  than  250  daily  observations  for  either  real 
or  simulated  trajectories  contain  no  information  and  are  in 
white. 


96 


0.15 


Integral  Time  Scale,  Zonal  Direction,  Drifter  Data 
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Figure  39.  Histograms  of  Lagrangian  integral  time  scales 
(days)  for  drifter,  "same  start",  and  "random  start" 
trajectories,  zonal  component.  The  time  domain  is  January  1, 
1993  to  December  31,  1997.  Bins  with  fewer  than  250  daily 
observations  for  either  real  or  simulated  trajectories  are 
not  represented. 
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Figure  40.  Histograms  of  Lagrangian 
(days)   for   drifter,   "same   start", 
trajectories,   meridional   component. 
January  1,  1993  to  December  31,  1997. 
250   daily   observations   for   either 
trajectories  are  not  represented. 
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2.    Integral  Length  Scales 

There  is  much  better  agreement  in  magnitude  between  the 
model's  integral  length  scales  and  those  of  the  drifters. 
The  ratio  of  model  length  scale  to  drifter  length  scale 
ranges  from  0.3  to  2.9  in  the  zonal  direction  and  from  0.3 
to  5.4  in  the  meridional  direction,  but  a  far  greater  number 
of  bins  show  values  close  to  unity  than  for  the  time  scales . 
The  ratio  plots  show  overestimation  of  zonal  length  scales 
occurring  mainly  south  of  30°N  and  north  of  50°N,  with  the 
intervening  band  showing  underestimation.  The  pattern  of 
ratios  in  the  meridional  direction  is  similar  to  that  of  the 
time  scales,  with  elevated  values  north  of  the  Gulf  Stream 
and  in  the  southeast  Subtropical  Gyre.  The  exaggeration  is 
more  pronounced  in  the  northern  zone,  while  underestimation 
is  most  evident  in  the  lower  western  basin. 

Again,   "random  start"   (not  shown)   and  "same  start" 

length   scales   show   similar   spatial   distributions   and 

magnitudes.  The  histograms  for  the  two  trajectory  schemes 

(figs.  47,  48)  show  a  good  degree  of  similarity,  and  are 

close  in  shape  to  those  of  the  drifters.  The  zonal  length 

scale   histograms   show   longer   tails   for   the   drifters, 

extending  to  165  km  compared  to  134  km  for  the  "same  start" 

simulation.  The  number  of  occurrences  below  2  0  km  is  about 

four  times  higher  for  the  drifters  than  for  the  "same  start" 

regime.  For  all  three  histograms,  most  length  scales  are 

contained  between  15  and  60  km.  The  degree  of  agreement  for 

the  meridional  histograms  is  even  better.   In  all  three 
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nearly  all  length  scales  fall  between  10  and  60  km;  the  tail 
behavior  is  almost  identical,  with  extreme  values  reaching 
about  100  km.  The  only  significant  differences  are  skews  in 
the  model  histograms  not  seen  in  the  drifter  plot,  and 
distribution  differences  in  the  range  50  -  65  km. 
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Figure  41.  Lagrangian  integral  length  scale  (km),  zonal 
component,  drifters.  The  time  domain  is  January  1,  1993  to 
December  31,  1997.  Bins  with  fewer  than  250  daily- 
observations  contain  no  information  and  are  in  white. 
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Figure  42 .  Lagrangian  integral  length  scale  (km) ,  zonal 
component,  "same  start"  simulation.  The  time  domain  is 
January  1,  1993  to  December  31,  1997.  Bins  with  fewer  than 
250  daily  observations  contain  no  information  and  are  in 
white . 
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Figure  43.  Lagrangian  integral  length  scale  (km),  meridional 
component,  drifters.  The  time  domain  is  January  1,  1993  to 
December  31,  1997.  Bins  with  fewer  than  250  daily 
observations  contain  no  information  and  are  in  white. 


103 


"1 


""V  - ,  r 


\ 

,;' 

> 

1 

\ 

V 

" 

V 

^vJ. 

(■> 


o 


o 

CD 


o 

LO 


o 
■3- 


o 

CO 


o 

CM 


Figure  44.  Lagrangian  integral  length  scale  (km) ,.  meridional 

component,  "same  start"   simulation.   The   time  domain   is 

January  1,  1993  to  December  31,  1997.  Bins  with  fewer  than 

2  50  daily  observations  contain  no  information  and  are  in 

white . 
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Figure  45.  Ratio  of  "same  start"  Lagrangian  integral  length 
scale  to  drifter  integral  length  scale,  zonal  direction.  The 
time  domain  is  January  1,  1993  to  December  31,  1997.  Bins 
with  fewer  than  250  daily  observations  for  either  real  or 
simulated  trajectories  contain  no  information  and  are  in 
white . 
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Figure  46.  Ratio  of  "same  start"  Lagrangian  integral  length 
scale  to  drifter  integral  length  scale,  meridional 
direction.  The  time  domain  is  January  1,  1993  to  December 
31,  1997.  Bins  with  fewer  than  250  daily  observations  for 
either  real  or  simulated  trajectories  contain  no  information 
and  are  in  white. 
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Figure  47 .  Histograms  of  Lagrangian  integral  length  scales 
(km)  for  drifter,  "same  start",  and  "random  start" 
trajectories,  zonal  component.  The  time  domain  is  January  1, 
1993  to  December  31,  1997.  Bins  with  fewer  than  250  daily 
observations  for  either  real  or  simulated  trajectories  are 
not  represented. 
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Figure  48.  Histograms  of 
(km)   for   drifter,   "same   start", 
trajectories,   meridional   component. 
January  1,  1993  to  December  31,  1997. 
250   daily   observations   for   either 
trajectories  are  not  represented. 
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3 .    Comparison  of  Integral  Time  and  Length  Scales  with 
Previous  Studies 

The  integral  time  scales  computed  from  the  drifter  data 

compare  well  with  a  previous  large-area  study  by  Krauss  and 

Boning   (19  87),   those   from   corresponding   areas   agreeing 

generally  within  1  day.  Integral  length  scales  show  less 

agreement.  The  exact  values  are  given  in  the  table  below: 


LONGITUDE 

LATITUDE 

ITSZ 
K&B 

ITSZ 
WOCE 

ITSM 
K&B 

ITSM 
WOCE 

ILSZ 
K&B 

ILSZ 
WOCE 

ILSM 
K&B 

ILSM 
WOCE 

3  0°N- 
33°N 

16°W-27°W 

4.5 

4.6 

5.6 

2.9 

32 

49 

40 

29 

39°N- 
43°N 

22.5°W- 
29°W 

8.6 

6.3 

5.7 

5.8 

80 

53 

53 

50 

45°N- 
47.5°N 

2  0°W- 
2  8.5°W 

6.0 

6.2 

3.4 

2.7 

72 

72 

41 

33 

3  9°N- 
43.5°N 

38.5°W- 
43 .5°W 

3.0 

3.0 

2.3 

1.5 

53 

85 

42 

40 

43°N- 
47°N 

3  2.5°W- 
41.5°W 

2.8 

2.4 

3.2 

1.7 

56 

54 

66 

35 

47°N- 
52.5°N 

31.5°W- 
43 .5°W 

2.1 

1.8 

2.1 

1.6 

45 

36 

47 

34 

48.5°N- 
52.5°N 

20°W- 
3  0.5°W 

3  .2 

4.5 

2.2 

2.5 

59 

63 

37 

33 

53°N- 
55°N 

3  0°W- 
37.5°W 

4.3 

4.3 

3.4 

2.4 

58 

36 

46 

21 

Table  1.  Comparison  of  Lagrangian  integral  time  and  length 
scales  with  those  of  Krauss  and  Boning  (1987)  . 
Abbreviations  used  are  ITSZ  (integral  time  scale,  zonal) , 
ITSM  (integral  time  scale,  meridional),  ILSZ  (integral 
length  scale,  zonal),  ILSM  (integral  length  scale, 
meridional) ,  K&B  (Krauss  and  Boning) ,  WOCE  (World  Ocean 
Circulation  Experiment  (with  which  the  drifter  data 
collection  is  associated) ) . 


Time  and  length  scales  are  found  by  Colin  de  Verdiere 

(1983)  in  the  Bay  of  Biscay  to  be  2.5  days  and  22.3  km  for 

the   zonal  direction,   and  2   days   and  17.1  km  for   the 

meridional  direction.  For  roughly  the  same  region  we  obtain 
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4.3  days  and  32.1  km  for  the  zonal  direction  and  2.6  days 
and  21.5  km  for  the  meridional  direction.  Again,  the  time 
scales  show  better  agreement  than  the  length  scales. 

4.    Eddy  Dif fusivities 

To  estimate  eddy  dif fusivities,  we  take  the  maximum 
values  attained  by  the  integrals  of  the  Lagrangian 
autocovariance  functions  (see  Appendix  A)  over  the  time 
domain  0-25  days.  This  is  a  reasonable  simplification,  and 
for  an  ensemble  of  particles  released  simultaneously  from 
one  point  into  a  perfectly  homogeneous,  stationary 
turbulence  field,  yields  a  value  similar  to  the  theoretical 
result,  given  that  the  integral  time  scale  is  generally  much 
smaller  than  25  days. 

The  limitations  of  the  method  should  be  made  clear, 
however.  The  assumption  of  homogeneity  is  unrealistic  in 
some  locations,  as  stated  previously.  Strong  inter-annual 
and  intra-annual  effects  in  the  velocity  field  violate  the 
stationarity  assumption.  Such  departures  from  the  ideal 
state  lead  to  autocovariance  functions  that  do  not  asymtote 
smoothly  to  zero;  the  resulting  eddy  dif fusivities  may  have 
significant  oscillations  that  complicate  specifying  a 
representative  value. 

Difficulties   in   quantifying   eddy   diffusivity   are 

exacerbated  by  large  time  lags.  Even  if  a  particle  begins 

its  path  in  a  locally  homogeneous  and  stationary  field,  its 

track,  on  which  the  diffusivity  calculation  is  based,  may 
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lead  into  a  region  of  significantly  different  properties 
over  the  time  lag  of  interest.  Further,  actual  computations 
of  eddy  diffusivity  based  on  field  data  cannot  meet  the 
requirement  of  simultaneous  particle  release  from  one  point. 
The  point  of  release  must  be  approximated  by  a  two- 
dimensional  region;  the  tracks  falling  into  that  area  are 
then  used  to  compute  an  average  Lagrangian  autocovariance 
function,  which  is  then  integrated  to  obtain  the 
diffusivity.  This  highlights  another  difficulty  with  large 
time  lags.  Because  drifter  tracks  are  finite,  the  number  of 
contributions  to  the  average  autocovariance  function 
decreases  with  increasing  time  lag  as  more  track 
terminations  occur. 

In  regions  where  velocities  have  high  temporal 
autocorrelation,  the  integrals  will  continue  to  increase  up 
to  large  time  lags.  This  forces  a  compromise:  one  must 
either  cite  a  diffusivity  that  clearly  underestimates  the 
steady-state  value,  or  else  incur  the  errors  discussed  above 
in  pursuing  higher  time  lags.  Examples  of  such  behavior  in 
the  drifter  data  include  the  Gulf  Stream,  for  which  there  is 
strong  zonal-zonal  autocorrelation  out  to  25  days,  and  north 
of  50°N  and  west  of  45°W,  where  both  diagonal  components 

( /^(T)  ,  /q2  (1)  )     of     diffusivity     indicate     persistent 

autocorrelation.   The  model  shows  the  same  tendencies  in 
these  areas;  additionally,   for  most  of  the  band  2  0°N  to 
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3  0N°,  k1x(t)      shows  no  sign  of  reaching  a  plateau  over  2  5 

days.  The  effect  is  most  pronounced  in  the  west  basin. 

The  above  discussion  illustrates  the  difficulty  in 
choosing  a  practical  rule  for  determining  eddy  diffusivity 
values.  While  our  taking  a  simple  maximum  value  over  a  fixed 
time  domain  may  seem  naive,  it  should  give  accurate  results 
in  "well-behaved"  regions.  Where  the  dif fusivities  fail  to 
plateau,  the  estimates  must  be  interpreted  as  lower  bounds. 
In  more  ill-behaved  areas  the  use  of  more  complex  methods 
based  on  asymptotic  diffusion  theory  seems  questionable, 
given  the  departure  from  either  homogeneity,  stationarity, 
or  both. 

Figures  49  and  50  display  the  integrals  of  the  sample 
Lagrangian  velocity  autocovariance  functions  from  four 
distinctive  regions  of  the  North  Atlantic,  for  the  drifter 
data  and  "same  start"  trajectories.  Only  zonal-zonal  and 
meridional-meridional  components  are  shown.  The 
characteristics  of  the  regions  are  as  follows: 

Region  1  (55°N-60°N,  35°W-40°W) :  This  region  is  located 
southeast  of  the  cusp  of  Greenland.  It  is  representative  of 
latitudes  north  of  55°N,  where  diffusivity  magnitudes  are 
generally  within  5xl07  cm2/ sec,  plateaus  are  reached  within 
25  days  in  most  cases,  inhomogeneities  are  not  pronounced, 
and  where  there  is  generally  good  agreement  between 
diffusivity  curves  from  the  drifters  and  from  the  "same 
starts"  trajectories. 

112 


2.5 


E 
u 

O 


>1.5 


0.5 


Drifter  Data,  Region  1 
55N-60N,  35W-40W 


O^r 


-j.Zi 

■       A 

- 

A 

■ 

■ 

A  =  zonal 

■+■  =  meridional 

. 

10  20 

Time  Lag  (Days) 


30 


'Same  Starts',  Region  1 
55N-60N,  35W-40W 


10  20 

Time  Lag  (Days) 


30 


25  r 

1=  20 
o 

'o 

5. 

1-15 


Q 
>.10 


5^ 
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Figure  49.  Integrals  of  zonal  and  meridional  components  of 
the  sample  Lagrangian  velocity  autocovariance  functions 
(calculated  per  Appendix  A)  for  regions  1  and  2.  The 
diagonal  components  of  eddy  diffusivity,  on  which  the 
Lagrangian  integral  time  and  length  scales  are  based 
(equations  13,14),  are  taken  as  the  maxima  of  the 
corresponding  integrals  over  25  days. 
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Figure  50.  Integrals  of  zonal  and  meridional  components  of 
the  sample  Lagrangian  velocity  autocovariance  functions 
(calculated  per  Appendix  A)  for  regions  3  and  4.  The 
diagonal  components  of  eddy  diffusivity,  on  which  the 
Lagrangian  integral  time  and  length  scales  are  based 
(equations  13,14),  are  taken  as  the  maxima  of  the 
corresponding  integrals  over  2  5  days. 
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Region  2  (40°N-45°N,  60°W-65°W)  :  This  box  is  to  the 
south  of  Nova  Scotia,  in  the  path  of  the  Gulf  Stream.  It 
characterizes  the  western  North  Atlantic  between  latitudes 
35°N  and  45°N.  This  entire  region  is  under  direct  Gulf 
Stream  influence.  Inhomogeneities  are  evident  in  most  of  the 
drifter  diffusivity  plots.  Zonal  diffusivity  magnitudes  are 
significantly  larger  than  meridional  values,  and  quite  large 
in  an  absolute  sense  (exceeding  5x108  cm2/sec  for  the 
drifters)  .  In  about  half  the  cases,  for  both  drifters  and 
model,  zonal  dif fusivities  do  not  plateau  within  25  days. 
Model  dif fusivities  are  smaller  than  those  of  the  drifters, 
and  the  zonal  magnification  is  absent  in  the  model  over  the 
southern  half  of  the  region. 

Region  3  (35°N-40°N,  15°W-20°W) :  This  is  an  area  midway 
between  the  Azores  and  the  southern  coast  of  Spain.  It 
characterizes  most  of  the  lower  east  between  latitudes  40 °N 
and  3  0°N,  and  most  of  the  west  basin  between  3  0°N  and  3  5°N. 
Both  diagonal  components  of  diffusivity  are  distinctly 
underestimated  in  this  region.  Drifter  dif fusivities  are 
moderate,  ranging  from  about  0.5xl07  cm2/s2  to  2.0xl08  cm2/s2. 
Both  model  and  drifter  dif fusivities  tend  to  plateau  within 
25  days.  Inhomogeneities  are  seen  in  the  drifter 
dif fusivities  at  the  extreme  west  and  east  boundaries, 
associated  with  the  Gulf  Stream  and  the  Canary  Current, 
respectively. 

Region  4  (35°N-40°N,  15°W-20°W) :  This  region  is  located 

at  the  edge  of  the  North  American  Basin,  paralleling  the 
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Florida  peninsula.  It  is  representative  of  gradually  rising 
dif fusivities  from  east  to  west  in  the  latitude  band  20°N- 
30°N,  but  is  distinct  from  the  Gulf  Stream  region  to  the 
north  in  that  the  two  components  of  diffusivity  are  roughly 
of  the  same  magnitude,  and  the  effects  of  inhomogeneities 
are  not  strong.  The  model  underestimates  dif fusivities, 
particularly  the  meridional  component,  in  this  region,  and 
excessive  zonal  autocorrelation  is  indicated  in  the  model  by 
the  absence  of  plateaus. 

As  indicated  by  equation  13 ,  the  relative  magnitudes  of 
the  eddy  dif fusivities  seen  in  Figures  49  and  50  are  related 
to  the  product  of  the  zero-lag  autocovariances  and  the 
Lagrangian  integral  time  scales.  The  zero-lag 
autocovariances  are  mathematically  equivalent  (and 
computationally  nearly  equivalent)  to  the  Eulerian 
variances . 

In  region  1  the  integral  time  scales  (see  Figures  37, 

38)  are  seen  to  be  overestimated  by  a  factor  of  about  two. 

Meridional  variances  appear  to  be  slightly  underestimated  by 

the  model   (Figures  18  and  19)  while  zonal  variances  are 

about  half  those  of  the  model.  This  is  consistent  with  the 

top  two  plots  of  Figure  49,  which  show  nearly  equal  zonal 

components   of   diffusivity,   while   the   model   slightly 

exaggerates  the  meridional  component.   In  region  2,   zonal 

integral  time  scales  are  overestimated  by  a  factor  of  about 

two,   while  meridional  scales  are  roughly  the  same.   Both 

components  of  variance  are  underestimated,  but  mainly  in  the 
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zonal  direction.  The  result  is  that  both  diffusivity 
components  are  slightly  smaller  in  the  model,  as  shown  in 
the  bottom  graphs  of  Figure  49.  Region  3  shows  model 
dif fusivities  to  be  significantly  low  (Figure  50).  Time 
scales  are  nearly  equal  between  model  and  data,  but  the 
distinctively  low  model  variances  in  the  eastern  basin  cause 
its  dif fusivities  to  be  underestimated  by  over  three  times, 
as  seen  in  the  top  plots  of  Figure  50.  Finally,  in  region  4 
the  zonal  time  scale  ratio  is  again  about  two  and  the 
meridional  ratio  nearly  unity.  Both  variances  are 
underestimated  roughly  by  a  factor  of  two,  giving  the 
expected  result  in  dif fusivities,  as  indicated  by  the  bottom 
plots  of  Figure  50. 

Though  not  shown,  eddy  diffusivity  plots  from  the 
"random  start"  trajectories  agree  closely  with  those  of  the 
"same  start"  scheme  over  the  entire  study  domain,  reflecting 
the  same  patterns  with  respect  to  the  drifter  dif fusivities . 
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VI .   CONCLUSIONS 

Both  Eulerian  and  Lagrangian  comparisons  are  used  in 
this  study  to  explore  the  adequacy  of  the  1/6-degree  POP 
model  in  representing  the  general  surface  circulation  of  the 
North  Atlantic  Ocean.  A  preliminary  Eulerian  analysis  of  the 
1/10-degree  model  is  done  for  a  two-year  subset  of  the 
period  of  the  original  data.  The  findings  are  in  general 
agreement  with  previous  studies  of  the  North  Atlantic 
circulation  and  comparisons  of  the  model  with  various  data 
sources . 

The  salient  shortcomings  of  the  1/6-degree  POP  model 
are  its  under-representation  of  flow  variability  and 
spurious  dynamic  effects  in  the  upper  Gulf  Stream  and 
Subtropical  Gyre.  These  in  combination  are  the  source  of  the 
discrepancies  seen  in  comparison  with  our  drifter  data  set. 

The  Eulerian  comparisons  show  the  1/6-degree  model  to 
resolve  the  major  flow  features  indicated  by  the  drifters 
over  the  North  Atlantic  Ocean.  The  northern  marginal 
currents  are  especially  well  modeled.  Exceptions  to  the 
general  agreement  include  the  model's  inaccurate 
representation  of  major  direction  changes  in  the  Gulf  Stream 
Extension,  unrealistically  high  structure  in  the  Subtropical 
Gyre's  recirculation  patterns,  poor  resolution  of  the  Canary 
Current,  and  apparent  inability  to  resolve  the  Azores  Front. 
The  effects  of  the  model's  artificial  front  in  the  southern 
Subtropical  Gyre  and  anticyclone  in  the  Gulf  Stream  between 
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Cape  Hatteras  and  Cape  Cod  are  clearly  evident  in  the 
Eulerian  plots. 

The  model's  total  contribution  to  MKE  appears  to 
be  correct;  drifter  and  model  MKE  histograms  show  almost 
identical  frequency  distributions  as  well.  The  spatial 
distribution  of  this  energy  seems  to  be  perturbed  by  the 
model's  artifacts  and  its  inaccurate  flow  geometries  in  the 
Gulf  Stream  Extension,  however.  Energy  is  excessive  where 
the  Gulf  Stream  closely  follows  the  continental  shelf,  and 
on  the  eastern  side  of  the  Mid  Atlantic  Ridge,  where  the 
model's  Northwest  Corner  occurs.  Elsewhere  in  the  Gulf 
Stream  MKE  is  underestimated.  There  is  both  overestimation 
and  underestimation  in  the  Subtropical  Gyre,  associated  with 
the  model's  spurious  front  and  its  unresolved  Azores 
Current,  respectively.  No  significant  tendency  of 
underestimation  or  overestimation  by  the  model  as  a  function 
of  energy  level  is  found. 

The  1/6-degree  model  appears  to  nearly  uniformly 

under-represent  variance  about  the  surface  mean  flow,  though 

the  singularities  of  the  model  both  amplify  and  diminish  the 

contrast.  The  model's  spurious  anticyclone  gives  elevated 

values   in   the   region   34°N-38°N,   70°W-76°W.   The   North 

Atlantic  Current  and  northeast  drift  north  of  the  Rockall 

Plateau  appear  closely  matched  in  EKE.   This   result   is 

misleading,  though,  because  the  model's  displacement  of  the 

Northwest    Corner    produces    reciprocal    effects    of 

underestimation    and    overestimation    west    and    east, 
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respectively,  of  the  Mid  Atlantic  Ridge.  Comparisons  in 
other  areas  are  not  complicated  by  the  model's  spurious 
effects;  the  continental  boundary  currents  in  the  northeast 
Atlantic  Ocean,  for  example,  are  unambiguously  well  modeled, 
while  EKEs  are  clearly  underestimated  over  most  of  the 
Subtropical  Gyre  and  in  the  central  Labrador  Sea. 

A  frequency  histogram  shows  a  large  overabundance  of 
model  EKE  values  below  50  cm2/s2  and  too  few  above  1500 
cmVs2.  As  with  MKE,  the  model  shows  no  strong  trend  toward 
underestimating  or  overestimating  EKE  as  a  function  of 
energy  level . 

The  model  represents  eddy  dif fusivities  fairly  closely 
at  latitudes  north  of  50N.  With  the  exception  of  the  Gulf 
Stream,  dif fusivities  are  distinctly  underestimated  by  the 
model  south  of  50°N.  Zonal  dif fusivities  are  also  markedly 
underestimated  in  the  Gulf  Stream.  These  patterns,  in 
conjunction  with  the  model's  nearly  uniformly  low  variances, 
determine  the  spatial  distribution  of  the  integral  length 
and  time  scales. 

Inhomogeneities  are  most  evident  along  the  course  of 

the  Gulf  Stream,   and  to  a  much  lesser  extent  the  other 

continental  boundary  currents;  they  generally  appear  less 

pronounced  in  the  model  than  in  the  drifter  data.  Failure  of 

dif fusivities  to  reach  plateaus  within  25  days,  indicating 

persistent  temporal  autocorrelation,  is  pronounced  for  the 

zonal-zonal  component  of  model  diffusivity  over  most  of  the 

mid  and  lower  east  basin  and  lower  west  basin. 
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The  model's  integral  time  scales  are  high,  by  average 
factors  of  about  1.9  and  1.8  in  the  zonal  and  meridional 
directions,  respectively.  Spatially,  the  overestimates 
occurred  mostly  in  northern  latitudes,  with  the  Gulf  Stream 
providing  a  demarcation  over  most  of  the  northern  basin. 
Zonal  time  scales  are  particularly  overestimated  in  the 
latitude  band  range  20°N  -  30°N.  Meridional  time  scales  are 
overestimated  by  factors  upwards  of  three  in  latitudes  south 
of  30°N. 

Length  scales  are  more  closely  matched;  they  are 
underestimated,  but  only  by  average  factors  of  1.1  and 
slightly  greater  than  unity  in  the  zonal  and  meridional 
directions,  respectively.  Again,  the  Gulf  Stream  divides 
higher  values  to  the  north  and  lower  values  to  the  south. 
The  zonal  length  scale  plots  show  the  same  band  of  elevated 
values  in  the  latitude  range  20°N  -  30°N.  A  significant 
exception  to  the  pattern  of  high  time  and  length  scales  to 
the  north  is  that  the  East  Greenland  Current  is  consistently 
well  represented  by  the  model. 

For  both  Lagrangian  integral  time  and  length  scales, 
spatial  patterns  shown  by  the  "random  start"  and  "same 
start"  regimes  are  similar;  frequency  histograms  also  tally 
fairly  closely  for  both  regimes.  We  evidently  do  not  incur 
substantial  bias  by  sampling  the  model's  velocity  fields  in 
a  deterministic  manner. 
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Preliminary  Eulerian  comparisons  of  the  NOGAPS  model 
run  with  the  drifter  data  over  the  period  January  1,  1993  - 
December  31,  1994  show  a  substantial  improvement  in  EKE 
representation.  The  turbulence  of  the  Subtropical  Gyre  is 
very  well  represented  in  the  NOGAPS  run,  where  the 
1/6-degree  model  shows  far  too  much  structure  in  the  surface 
flow.  The  Azores  Current  appears  in  the  1/10-degree  model 
output,  a  result  also  noted  by  Smith  et  al .  (1999),  and  the 
spurious  front  at  latitudes  2  6-N  -  3  02N  is  much  attenuated. 
The  Northwest  Corner  appears  in  the  correct  longitude  band, 
and  the  model's  artificial  anticyclone  seems  to  be  absent. 
There  continues  to  be  little  evidence  for  the  Canary  Current 
in  the  1/10-degree  run,  but  this  may  be  the  result  of  the 
sparsity  of  the  data;  the  current  is  also  more  poorly 
represented  in  the  two-year  subset  of  drifter  data  than  in 
the  full  drifter  data  set. 

Difficulties   remain   in   connection  with   the   energy 

balance   for   the   NOGAPS   formulation.   While   the   overall 

distribution  of  EKE  shows  a  large  improvement  over  that  of 

the  1/6-degree  model,  the  different  formulation  apparently 

has  the  attendant  effect  of  exaggerating  MKE .  Energies  above 

500  cm2/s2  are  much  more  frequent  than  in  the  1/6-degree 

model.  Though  the  large  surplus  of  EKE  in  the  0-50  cmVs2 

range   is   corrected   in   the   1/10-degree   model,    an 

overabundance   of   energy   appears   above   1000   cm  /s  . 

Geographically,  there  appears  to  be  a  meridional  mismatch  in 

EKE  along  the  course  of  the  Gulf  Stream.  The  gradient  of 
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EKE,  in  the  "stream  coordinate"  sense,  is  too  large  along 
the  southwest  boundaries  of  the  North  American  Basin.  The 
region  of  elevated  EKE  where  the  Gulf  Stream  breaks  from  the 
continental  shelf  extends  too  far  to  the  south  and  east; 
likewise,  the  region  of  high  EKE  associated  with  the 
Northwest  Corner  has  an  excessive  peripheral  reach. 

The  low  variability  of  surface  flow  in  the  1/6-degree 
model  is  probably  attributable  to  its  relatively  low  grid 
resolution  and  lack  of  a  surface  mixed  layer.  This  study 
makes  no  attempt  to  discern  the  relative  degrees  of 
improvement  in  adding  a  mixed  layer  and  increasing 
resolution.  However,  comparisons  of  sea  level  height  data 
with  a  previous  1/10-degree  run  using  ECMWF  winds  and  no 
mixed  layer  (Smith  et  al .  ,  1999)  show  improvement  on  about 
the  same  order  as  that  which  we  observe.  Resolution  is 
therefore  probably  the  dominant  factor. 

Evidently  more  experimentation  is  needed  in  order 
to  achieve  a  more  accurate  distribution  of  EKE  in  the 
1/10-degree  model,  particularly  to  attenuate  the  higher 
energies  while  still  elevating  energies  out  of  the  0  to  50 
cmVs2  range.  Ideally,  this  would  be  done  in  a  way  that  does 
not  inflate  MKE  on  the  whole,  as  apparently  occurs  in  the 
present  version.  Varying  forcing  frequencies,  and  friction 
coefficients  may  achieve  an  adequate  decoupling  of  the  two 
effects . 

On  the  basis  of  Eulerian  statistics,   the  1/6-degree 

model  appears  to  perform  fairly  well  at  latitudes  above 
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60°N.  Northern  boundary  currents  are  particularly  well 
modeled.  The  major  flow  structures  of  the  North  Atlantic  are 
well  described  by  the  1/6-degree  model  over  much  of  the 
domain,  with  the  exceptions  of  the  Gulf  Stream  Extension, 
the  spurious  anticyclone,  and  the  Subtropical  Gyre 
recirculation.  Outside  of  these  regions,  the  model  mainly 
suffers  from  generally  low  variability,  but  even  this  is  not 
pronounced  over  much  of  the  upper  latitudes.  The  Lagrangian 
integral  time  and  length  scales,  on  the  other  hand,  tend  to 
be  overestimated  in  the  upper  latitudes,  except  for  the 
southwest  continental  margin  of  Greenland. 

Our  overall  assessment  is  that  the  generally  low 
variability,  limited  feature  resolution,  displaced  branching 
of  the  North  Atlantic  Current,  and  spurious  circulation 
effects  observed  in  the  1/6-degree  model  make  it  inadequate 
for  use  in  an  integrated  global  forecast  system.  Based  on 
our  preliminary  results  from  the  1/10-degree  run,  we  find 
this  higher-resolution  version  to  substantially  correct  the 
problems  observed  in  the  1/6-degree  model.  It  appears  that 
the  1/10-degree  model  may  well  satisfy  the  requirements  of 
the  coupled  global  forecast  system,  and  accordingly  we 
recommend  further  testing  at  this  resolution. 
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APPENDIX  A.  CALCULATION  OF  THE  SAMPLE 
LAGRANGIAN  VELOCITY  AUTOCOVARIANCE  FUNCTION 

The  following  approach  is  that  of  Poulain  et  al . 
(1996).  As  described  in  section  V.A,  the  Lagrangian  mean 
velocity  of  an  ensemble  of  particles  is  referenced  to  their 
deployment  from  a  single  point  x  at  a  given  time  t0 .  In 
practice,  buoys  are  not  deployed  from  one  location  but 
rather  across  an  extended  geographical  area.  Therefore,  to 
obtain  the  sample  Lagrangian  velocity  autocovariance  for  a 
collection  of  actual  buoys  in  a  given  neighborhood,  the 
spatial  domain  of  the  definition  is  expanded  (from  zero)  to 
2  dimensions.  Any  buoy  that  inhabits  the  2-dimensional  area 
at  any  time  during  its  entire  record  is  considered  to 
contribute  to  the  "ensemble"  of  deployments  within  the 
region.  There  is  one  "deployment"  corresponding  to  each 
time  step  spent  in  the  region,  however,  each  deployment  is 
considered  to  take  place  at  the  same  reference  time,  t0 
(hence,  we  assume  stationarity) .  For  our  purposes,  the 
specific  value  of  t0  is  irrelevant  and  set  to  zero. 

Consider  a  buoy  whose  displacements  from  an  arbitrary 
point  are  measured  over  a  discrete  time  interval  [m,n]  as 
(rm,rm+1,rm+2,...,rn)  .  Let  some  contiguous  segment  ( r. ,  ri+1 ,  ri+2 , ... ,  r . ) 
of  its  track  pass  through  the  region  of  interest,  and 
suppose  the  sample  Lagrangian  autocovariance  function  over  a 
time  lag  interval   [-T,T]   is  sought  for  that  region.  Then 
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that  segment's  contribution  to  the  function  is  derived  from 
the  series  of  overlapping  "pseudotracks"  given  by 
(rr    r     rrrrr...r    r    ri 

\  i-r'   i-if-li'  i-ir-2) i-2'  i-V  i'  i+V  i+2 i+(T-2)'  i+<r-l)'   i+T/ 

rr  r  rrrrr...r  r  ri 

\  [i+l)-T     (i+l)-(T-l)     (i+l)-<T-2)         '    ii+h-2'   (i+l)-l     i+1     (i+l)+l     (i+l)+2        '   (i+l)+(T-2)     (i+l)+(T-l)     (i+h+r  " 

rr  r  rrrrr...r  r  r 

V   (J-l)-X'     (J-l)-(T-l)'     (J-l)-(T-2)'  '     (J-l)-2       (i-l)-l'     J-1'     (j-l)+]       (j-l>+2  (J-l)+(T-2)       (j-l)+(T-l)       (J-l)+T/' 

rr         r  rrrr  r         r         r 

V  j-T'   j-i.T-1)'  j-(T-2) j-2'   i-l'j"  j+V   j+2 j+.(T-2}'   j+(T-l)'  j+r/ 

where  velocity  information  is  available  for  all  pseudotrack 
elements  whose  time  indices  are  within  the  time  domain  [m,n] 
of  the  track.  Velocities  for  elements  whose  indices  fall 
outside  the  time  domain  are  assigned  null  values.  This 
process  is  repeated  for  all  track  segments  that  pass  through 
the  region,  generating  one  velocity  time  series  for  each 
observation.   For  the  kth  discrete  lag  xk  in   [-T,x],   the 

"ensemble   average"   V(Tk,t0,x)   =   <v(  t0+Tk  |x,  t0)  >L   is   then 

computed   as    the   mean   of    the   non-null   velocities 
corresponding  to  the  kth  pseudotrack  elements.  The  residuals 

v'  ( t0+xjx,  t0)   about  the  kth  lag  averages  follow  directly, 

where  the  zero-lag  residual  is  the  special  case  v'(t0|x,t0). 
The  value  of  the  sample  Lagrangian  autocovariance  function 

at  time  lag  Tk  for  components  i  and  j  (where  i  and  j  may  be 

either  zonal  or  meridional)  is  then  given  by  the  mean  of  the 

cross-products  vi '  ( t0+Xk  |x,  t0)  v/(t0|x,t0). 
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APPENDIX  B.  NUMERICAL  GENERATION  OF  PARTICLE  TRAJECTORIES 
USING  THE  RUNGE-KUTTA  METHOD 


The  generation  of  numerical  particle  trajectories  for 
the  "random  start"  and  "same  start"  schemes  involves  a 
series  of  initial  value  problems  to  the  first  order 
differential  equation 

v(t)=/(r,0,  r(0)  =  r0. 

Here,  r  is  a  displacement  vector  defined  over  the  spatial 
domain  of  the  model,  and  v  is  a  velocity  vector  obtained  by 
bilinear  interpolation  of  the  model  velocity  field  at  r. 
Given  displacement  r(t)  at  time  t,  the  displacement 
r(t+At)  at  time  t  +  At,  is  given  by 


;+A/ 


r(t  +  At)  =  r(t)+   Jv(r,/V'- 


For  the  model's  discrete  time  step  At,   a  practical 

algorithm  for  advancing  a  solution  from  r(t)   to  r(t+At) 

employs  the  fourth  order  Runge-Kutta  method,  described  in 
detail  in  Press  et  al .  (1992).  The  method  derives  from  the 
simpler  Euler  Method,  which  uses  the  first  term  of  the 
Taylor  expansion  of  r  about  t.  In  the  context  of  our  problem 
this  results  in  the  familiar  formula 

r(t  +  At)  =  r(t)  +  v(r,  t)At. 
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The  disadvantage  of  the  Euler  method  is  that  it  uses 
derivative  information  (v)  only  at  the  start  of  the  time 
interval.  The  method  has  an  error  only  one  power  smaller 

than  the  error  of  the  time  step,  or  0((At)2)  .  A  more  stable 
and  accurate  method  uses  the  derivative  at  time  t  to  take  a 
step   to   the   midpoint   r(t+At/2)   of   the   interval.   The 

derivative  v(r(t+At/2),  t+At/2)  is  then  substituted  for  the 

original   v(r,t)    to   extrapolate   to   the   next   point; 
mathematically, 


k,  =v(r,r)Af, 

k?  =  v(r  +— k,,r  +— At), 
2  2 

r(t  +  At)  =  r(t)  +  k ,  +  0((  At)3 ). 


The  error  term  0((At)2)   is   eliminated  by  symmetry, 

giving  higher  accuracy.  The  fourth  order  Runge-Kutta  method 
(a  method  is  called  nth  order  if  its  error  term  is  of  order 
n+1)  uses  similar  symmetric  canceling  with  derivatives 
evaluated  at  the  initial  point,  twice  at  trial  midpoints, 
and  once  at  a  trial  endpoint.  The  result  is 
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k,  =\(r,t)At, 

k,  =  v(r  +— k.,f  +  — Af)Ar, 
2     2 

k3  =  v(r  +— k2,r  +— Af)Ar, 
k4  =  v(r  +  k3,r  +  Af)A?, 

r(r  +  A?)  =  r(0  +  — +  —  +  —  +  — +  0((A/)5). 
6   3   3   6 


The  method  is  illustrated  graphically  in  Figures  51  and  52 
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Figure  51.  The  Euler  and  Midpoint  methods  for  numerical 
integration  of  ordinary  differential  equations.  The  notation 
is  adapted  to  the  problem  of  generating  trajectories.  The 
dependent  variable,  r,  is  actually  two-dimensional.  Plot  (a) 
illustrates  the  Euler  method,  which  uses  derivative  (v(t)) 
information  only  at  the  lower  endpoint  of  the  time  interval 
to  advance  the  solution  to  the  next  point.  Plot  (b)  shows 
the  more  accurate  midpoint  method,  in  which  the  derivative 
at  the  lower  endpoint  is  used  to  take  a  trial  step  to  the 
midpoint  of  the  interval.  The  derivative  at  the  midpoint  is 
then  used  to  extrapolate  the  function  value  over  the 
interval.  After  Press  et  al.,  (1992). 
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Figure  52  .  The  Fourth  Order  Runge-Kutta  method,  used  to 
generate  simulated  trajectories  for  the  "same  start"  and 
"random  start"  schemes.  The  method  is  based  on  the  Euler  and 
Midpoint  methods,  but  uses  derivative  information  once  at 
the  lower  endpoint  of  the  interval,  twice  at  trial 
midpoints,  and  once  at  the  endpoint.  The  four  derivatives 
are  then  used  to  advance  the  solution  across  the  interval. 
After  Press  et  al .  (1992). 
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